Gravitational waves from BH-NS binaries: Effective Fisher matrices and parameter estimation using higher harmonics

Gravitational waves from BH-NS binaries: Effective Fisher matrices and parameter estimation using higher harmonics

Hee-Suk Cho    Evan Ochsner    Richard O’Shaughnessy    Chunglee Kim    Chang-Hwan Lee Department of Physics, Pusan National University, Busan 609-735, Korea Center for Gravitation and Cosmology, University of Wisconsin-Milwaukee, Milwaukee, WI 53211, USA Department of Physics, West Virginia University, PO Box 6315, Morgantown, WV 26505, USA
July 15, 2019
Abstract

Inspiralling black hole-neutron star (BH-NS) binaries emit a complicated gravitational wave signature, produced by multiple harmonics sourced by their strong local gravitational field and further modulated by the orbital plane’s precession. Some features of this complex signal are easily accessible to ground-based interferometers (e.g., the rate of change of frequency); others less so (e.g., the polarization content); and others unavailable (e.g., features of the signal out of band). For this reason, an ambiguity function (a diagnostic of dissimilarity) between two such signals varies on many parameter scales and ranges. In this paper, we present a method for computing an approximate, effective Fisher matrix from variations in the ambiguity function on physically pertinent scales which depend on the relevant signal to noise ratio. As a concrete example, we explore how higher harmonics improve parameter measurement accuracy. As previous studies suggest, for our fiducial BH-NS binaries and for plausible signal amplitudes, we see that higher harmonics at best marginally improve our ability to measure parameters. For non-precessing binaries, these Fisher matrices separate into intrinsic (mass, spin) and extrinsic (geometrical) parameters; higher harmonics principally improve our knowledge about the line of sight. For the precessing binaries, the extra information provided by higher harmonics is distributed across several parameters. We provide concrete estimates for measurement accuracy, using coordinates adapted to the precession cone in the detector’s sensitive band.

pacs:
04.30.–w, 04.80.Nn, 95.55.Ym

I introduction

Ground based gravitational wave detector networks (notably LIGO Abb04 () and Virgo virgo ()) are analyzing results of design-sensitivity searches for the signals expected from the inspiral and merger of double compact binaries. Abb08 (); Abb09 (). For the lowest-mass compact binaries , the response of the detector to a binary merger with arbitrary masses, spins, and even eccentricity is well understood, particularly given the detectors’ limited and low sensitive frequency band Buo03 (); Buo04 (); Dam04 (); Pan04 (); Kon05 (); Buo05 (); Kon06 (); Tes07 (); Han08 (); Hin08 (); Aru09 (); Buo09 (). Though this complicated signal encodes all information about the binary’s spacetime Rya95 (), the amount of accessible information depends on signal strength (or signal-to-noise ratio) Jar12 (). Strong signals permit high-precision tests of general relativity; fainter signals allow high-precision constraints on some binary parameters; while very faint, short signals may only constrain the binary’s mass. Qualitatively speaking, we can distinguish two configurations if they are separated by contours , where (defined in Sec. III) is the (normalized) ambiguity function or “overlap” and is the signal-to-noise ratio (SNR).

As higher-order corrections and new physics are added to our models for gravitational-wave signals, the functional dependence on various parameters (such as masses, spins and orientation angles) in the model grows in complexity. On scales , for astrophysically plausible , the ambiguity function is generally smooth. However, it may have more complicated fine-scale structure which may not be detectable for expected signal strengths. The Fisher matrix approach to estimating parameter errors is based on differentiating a waveform with respect to its parameters. These derivatives are defined in an infinitesimal patch of parameter space, and are thus measuring fine-scale structure, which could potentially be misleading about larger-scale, observable trends. This point is illustrated by Fig. 2, where we fit a quadratic through the ambiguity function. The ambiguity function changes shape, so fits to small () and large () regions of parameter space would give rather different estimates of posterior widths and thus parameter accuracy. Similarly, if the standard Fisher matrix were calculated via finite difference, step sizes on these scales would give different results, with the smaller step size giving a misleading estimate of parameter accuracy for a signal of expected strength.

In this paper we propose a simple effective procedure to identify relevant scales and parameter correlations, construct suitable “effective Fisher matrices”, and estimate ambiguity functions at low but nontrivial signal to noise ratio. To demonstrate this technique, we examine the signal from selected black hole-neutron star (BH-NS) binaries, described in Section II. In this paper we use all available knowledge about the (post Newtonian) waveform, adopting a complete model for the adiabatic quasicircular inspiral of precessing BH-NS binaries. In particular, we employ all available harmonics and amplitude corrections, introducing small but non-negligible changes to the ambiguity function. In Section III we introduce our unconventional effective approach to the local ambiguity function. Using those tools, in Section IV we construct and approximate the overlap for signals similar to each reference binary. Motivated by parameter estimation, we provide explicit expressions for the Fisher matrix, correlations, and marginalized uncertainties for each configuration. For our fiducial configurations, higher harmonics principally allow us to improve our knowledge of the binary orientation, providing fairly little additional information about intrinsic parameters for the amplitude scales of immediate astrophysical interest.

Our results are complicated by coordinate-dependent effects, notably extreme sensitivity to the reference frequency at which parameters are specified. We show the choice of reference frequency can reduce (or introduce) fine-scale structure into the ambiguity function, similar to the effect of higher harmonics. Our effective approach can partially compensate for ill-chosen coordinates, such as the coalescence phase or initial spins. To reduce but not completely eliminate these systematic effects, we express our results using parameters specified near the peak sensitivity of the detector (here, 100 Hz).

i.1 Context and prior work

Several studies of gravitational wave detection from merging binaries have employed amplitude-corrected waveforms and higher harmonics. Investigations of space-based interferometers, such as the Laser Interferometer Space Antenna (LISA), have historically used complete signal models, accounting for both spin and precession Por08 (). As higher harmonics have a small effect, however, most previous studies of ground-based interferometers have omitted them, emphasizing spin. When included, higher harmonics were explored alone for non-precessing signals. Higher harmonics can allow detection of signals otherwise inaccessible due to the detector’s limited bandwidth Aru07 (); Bro07 (). The relative amplitudes of higher harmonics can probe astrophysical mechanisms for generating non-circularity Key10 (). Finally, higher harmonics (and precession) are well-known to break degeneracies and improve sky localization, particularly for LISA Lan06 (); Kle09 (); Lan11 ().

Several authors have explored the local ambiguity function “beyond the Fisher matrix”, including higher-order correlation functions Val08 (); Vit10 () and projection effects due to the local shape of the signal manifold Val11 (). These methods still use explicit derivatives of the ambiguity function (via explicit derivatives of the signal) to construct their series approximations.

The Fisher matrix is often nearly or exactly singular, making inversion numerically challenging. Several authors have pointed out that a singular value implies an unconstrained parameter, limited only by the prior; see, e.g., Vallisneri Val08 (). In many cases, including those singular values addressed in the text, the singular value corresponds to a bounded parameter (e.g., an angle). The singular value simply indicates that parameter cannot be measured. In the phenomenological limit described in this paper, precisely zero eigenvalues never occur, unless a parameter is constrained by symmetry.

Our goal in this work is to understand the typical shape of the posterior for a noise realization and coordinates in the signal space. Using one notion of “typical” would produce an average posterior over all noise realizations. Such an average posterior, however, could be slightly wider than the posterior from any given noise realization. Instead, in this work we attempt to characterize the typical shape of any one noise realization. To do so, in effect we “transport” each posterior so their peak likelihoods lie at the same point in parameter space. In practice, our procedure amounts to ignoring noise-realization-dependent changes to the posterior.

Ii Simulations

ii.1 Amplitude-corrected precessing waveform

In this paper we construct the post Newtonian (pN) gravitational wave signal from a BH-NS binary using the lalsimulation SpinTaylorT4 code lal (), which is an implementation based on the waveforms described in Buo03 (); Buo04 (). This time-domain code solves for the orbital dynamics of an adiabatic, quasicircular inspiralling binary by using the so-called TaylorT4 method (see Buo09 () for an explanation of this and similar methods) of evolving the orbital phase and frequency supplemented with precession equations to track the motion of the spins and orbital plane Kid95 (). The orbital phase and frequency evolution includes non-spinning corrections to 3.5pN order and spin corrections to 2pN order. The precession equations are given to 2pN order. This binary evolution is terminated prior to merger, either when it reaches the “minimum energy circular orbit”, or when the orbital frequency ceases to increase monotonically.

At each time, the values of the gravitational wave polarizations measured by a distant observer can be constructed from the orbital phase, orbital frequency and the orientations of the spins and orbital plane. We can construct either the commonly used “restricted” (i.e. leading-order) polarizations which contain only the dominant second harmonic of the orbital phase, or we can construct amplitude-corrected polarizations which contain terms that oscillate at other harmonics of the orbital phase (and also higher-order corrections to the second harmonic). Expressions for the polarizations valid for quasi-circular, precessing binaries are currently known to 1.5pN order Aru09 (); Kid95 (); Wil96 (). Throughout this work, when we refer to amplitude-corrected waveforms, we mean that we use the 1.5pN accurate polarizations.

ii.2 Simulation coordinates

For the precessing binaries, LIGO-scale studies have been complicated by poor choice of coordinates, associated with the start of the waveform. The waveform generation code of the standard LIGO algorithm defines all the geometrical parameter values at the initial frequency (40 Hz for the initial LIGO and 10 Hz for the advanced LIGO), and evolve the binary system to get the full waveforms. Specifically, the orbit is described at some point, including the spin, orbital angular momentum vector, and the orbital phase. By contrast, the detector is more sensitive to higher frequencies. Allowing for the decreasing signal strength with frequency Bro12 (), the detector is most sensitive to the instantaneous binary configuration at 100 Hz (initial LIGO) and 40 Hz (advanced LIGO, e.g., see Fig. 2 in Bro12 ()). Motivated by this fact, we choose the reference frequency (), at which the instantaneous orientations of the spins and orbital plane are defined, to be 100 Hz.

One significant effect introduced by setting the reference frequency at 100 Hz is related to the orbital phase. The ambiguity function is dependent on how we choose the reference frequency. In appendix B, we describe an example and illustrate the significant effects in detail; see also Fig. 2.

Furthermore, following Brown et al.Bro12 (), we define the geometrical parameters to be angles between the radiation vector , total angular momentum axis and orbital axis as in Fig. 1. For comparison, other conventions specify some parameters using the line of sight, a vector pointing from our detector to the binary.

For non-spinning or aligned-spin binaries, the total angular momentum is parallel to the orbital angular momentum and the orbital axis is fixed. In effect, the conventional radiation frame is equivalent to the geometrical frame.

ii.3 Fiducial simulations and local coordinates

In the case of a non-spinning binary, the binary is specified by 9 parameters. In this work, we choose masses () as intrinsic parameters, the polar (inclination) and azimuthal (polarization) angles () of the orbital axis with respect to the radiation vector and an orbital phase as extrinsic parameters. Because we maximize the ambiguity function over the polarization, we need not take this parameter into explicit account henceforth. Remaining parameters are, the distance to the detector, sky position (two angles), and the coalescence time. The fiducial values of parameters are summarized in Table 1. Mass components () can be expressed by the symmetric mass ratio and chirp mass , we adopt these parameters in this work.

If the NS spin is assumed to be 0, the aligned-spin binary is specified by 10 parameters. 9 parameters are the same as the non-spinning case and the additional intrinsic parameter is dimensionless BH spin parameter . The fiducial values are also summarized in Table 1.

The waveform of the precessing binary can be defined by 12 parameters if the NS spin is assumed to be 0. In this work, we consider , , BH spin , and the opening angle of the precessing cone as intrinsic parameters, , , , and the orbital phase as extrinsic parameters. Because we maximize over the polarization angle , the parameter is eliminated from further consideration. Remaining parameters are the distance, sky position (two angles), the coalescence time. Throughout this paper the units are solar masses (for ); radians (for angles); or the natural dimensionless units (for ).

Motivated by Bro12 (), we adopt a challenging reference configuration, where the polarization along the line of sight oscillates between circularly polarized ( along ) and linearly polarized ( perpendicular to ). Furthermore, to explore the extent to which higher-order harmonics allow measurement of parameters that only weakly impact the signal, we consider two possible sets of initial conditions for along its precession cone. The fiducial values of the parameters are summarized in Table 2. For case1, the orbital axis is perpendicular to the radiation vector at 100 Hz, for case2 it is parallel to the radiation vector at 100 Hz. All the parameter values are the same between both cases except for .

ii.4 Fiducial network

We assume two identical interferometers placed perpendicular to the incident signal, which is the optimal sky position of the source. We also assume the two interferometers are oriented by related to one another, giving comparable sensitivity to both polarizations. For the incident waveforms, we assume a zero noise limit to understand how similar the signals are. While not realistic, they avoid introducing complexity of the signal due to the source sky position.

Iii distinguishing simulations

iii.1 Ambiguity function

In this work, we reorganize the two projections of the strain tensor and into a complex function:

 h(t)≡h+(t)+ih×(t). (1)

We coherently compare a fiducial signal , where indicates a fiducial source parameter set, to a nearby signal , with parameters , by a complex overlap Osh12 ()

 ⟨h0|h⟩≡2∫∞−∞dfSn(f)[~h0(f)~h(f)∗], (2)

where is the Fourier transform of and is a detector strain noise power spectrum. For simplicity, we adopt a semianalytic initial LIGO sensitivity Dam01 (); Aji09 (). As pointed out by Osh12 (), this complex-valued expression characterizes the ability of a network to distinguish signals. The real part of the complex overlap corresponds to a linear sum of the conventional real overlaps of the two gravitational wave polarizations:

 Re⟨h0|h⟩=(h0,+|h+)+(h0,×|h×), (3)

where indicates the conventional overlap of two real functions defined by

 (h0|h)≡4∫∞0dfSn(f)Re[~h0(f)~h(f)∗]. (4)

In appendix A, we summarize the differences between the real and complex overlaps.

We note that a change of the polarization angle, , simple causes a rotation of the argument of the complex wave strain function, . Thus it is trivial to find the value of which makes the complex overlap purely real, so that

 Im⟨h0|h′⟩ = 0, (5) Re⟨h0|h′⟩ = (h0,+|h′+)+(h0,×|h′×),

and the value of the complex overlap maximized over polarization angle is simply

 maxψ Re⟨h0|h⟩=|⟨h0|h⟩|. (6)

The complex overlap (like the real-valued overlap) can also be maximized over the coalescence time via an inverse Fourier transform as described in All05 (). In particular, one uses the fact that

 ~h(tc=t)=~h(tc=0)e−2πift (7)

and notes that the inverse Fourier transform of the complex overlap integrand in Eq. (2) will compute the complex overlap for all possible coalescence times of at once

 ⟨h0|h(tc=t)⟩≡2∫∞−∞dfSn(f)[~h0(f)~h(f)∗]e2πift. (8)

The (normalized) ambiguity function between two waveforms and is then defined as the complex overlap maximized over polarization angle and coalescence time,

 P(λ0,λ)≡maxtc,ψ|⟨h0|h⟩|√⟨h0|h0⟩⟨h|h⟩. (9)

Unless otherwise noted, all overlaps are maximized in time and polarization. This is different from maximizing over orbital phase ; see appendix A and B.

iii.2 Likelihood

The detector noise is assumed to be a stationary and Gaussian process. Given the detector output representing a real-valued signal in real-valued noise, the probability for the noise to have some realization is Fin92 ()

 p(N=N0)∝e−(N0|N0)/2. (10)

The posterior probability that the gravitational wave signal is characterized by the parameters , can be expressed by , where is the prior probability that the signal is characterized by , is the likelihood, which can be written by Fin92 ()

 L(S|λ)=C×exp[−12(S−H(λ)|S−H(λ))], (11)

where is a proportional factor which, for simplicity, we assume to be 1 in this work.

Since we consider the complex strain, by choosing the appropriate polarization angle we shall write the detector output for the detector 1 and 2.

 S1=H++N1,   S2=H×+N2, (12)

also

 s=S1+iS2,   h0=H++iH×,   n0=N1+iN2. (13)

The probability for the noise to have both realizations and is

 p(N1,N2) ∝ e−⟨N1|N1⟩/2e−⟨N2|N2⟩/2 = e−Re⟨N1+iN2|N1+iN2⟩/2.

Finally, using Eqs. (12 - III.2), Eq. (11) can be expressed by the complex signals:

 L(s|λ)=exp[−12Re⟨s−h(λ)|s−h(λ)⟩]. (15)

Substituting into this equation Cut07 (), the likelihood is

 L = exp[−12Re⟨n0+h0−h|n0+h0−h⟩] = exp[−12Re{⟨h0−h|h0−h⟩+2⟨n0|h0−h⟩ +⟨n0|n0⟩}],

where the second and third terms in the a square bracket depend on the noise realization. They shift the position of the maximum likelihood but only weakly change the shape of the likelihood curve. In the limit of high SNR, these noise-dependent terms can be neglected, so,

 L = exp[−12Re⟨h0−h|h0−h⟩] = exp[−12{⟨h0|h0⟩+⟨h|h⟩−2Re⟨h0|h⟩}].

This equation corresponds to the case where two detectors are placed to have the maximum response to the incident two polarizations [for the detector placement, see Section II D]. While, Eq. (11) corresponds to one detector response to one polarization; see appendix A.

Using Eqs. (6) and (9), and the SNR defined by , the log likelihood can be expressed by our complex overlap convention:

 lnL=−ρ2(1−P) (18)

where we assume the same strength for both signals.

For a given log likelihood, the scale of interest of the ambiguity function depends on the signal strength :

 1−P≤1ρ2. (19)

By approximately identifying the surface of the likelihood, this condition allows us to estimate the set of parameters which cannot be distinguished from with a signal amplitude of using a signal model and noise curve that produces an overlap .111More properly, the probability of having likelihood lets us create a confidence volume for any target confidence level. Because the probability depends sensitively on , we anticipate the edge of this confidence interval will depend weakly (e.g., as ) on the precise probability used to define the threshold.

iii.3 Fisher matrix

If is close to , we can write to the first order in the error

 h0−h∼∂h∂λiΔλi. (20)

So, in the limit of high SNR, the likelihood [Eq. (III.2)] is given as , where is

 Γij=Re⟨∂h∂λi∣∣∣∂h∂λj⟩. (21)

This definition is analogous to the standard Fisher matrix, except that it is derived from the complex overlap and therefore contains information about both polarizations. If we assume that the prior is uniform, the parameter estimation errors (i. e., the posterior probability density function) can be expressed by the Gaussian distribution

 p(Δλi)=Ne−ΓijΔλiΔλj/2 (22)

where is the corresponding normalization factor.

Using another expression relating the Fisher matrix to the log likelihood Jar94 (); Val08 () and Eq. (18), we can write

 Γij = −∂2lnL(λ)∂λi∂λj∣∣∣λi=λj=λ0 = ρ2∂2(1−P)∂λi∂λj∣∣∣λi=λj=λ0=ρ2^Γij,

where is the fiducial value of source parameter. Here, we define the normalized Fisher matrix .

For Gaussian noise and high SNR, the inverse of the Fisher matrix is the covariance matrix () of parameter errors. The measurement error () of each parameter and correlation coefficient () between two parameters are defined as

 σi=√Σii,   cij=Σij√ΣiiΣjj, (24)

The correlation coefficients are -independent but often sensitive to small changes in . Conversely, the measurement error is inversely proportional to . For the purposes of illustration, we adopt whenever we calculate .

iii.4 Relevant scales and effective approach

The Fisher matrix formally involves derivatives, i.e., infinitesimal variations of a parameter . In this work, we compute an effective Fisher matrix by considering finite variations on scales which give physically observable changes to the ambiguity function, .

To understand the variability on multiple scales we plot a one-dimensional ambiguity function of for the leading-order amplitude, non-spinning binary in Fig. 2. In this figure, the ambiguity function is calculated via Eq. (9), changing only and fixing all other parameters to be the same for both signals. For comparison, we plot quadratic fits222Since the posterior function is a normal distribution in the limit of high SNR (see Eq. (22)), a quadratic fitting function usually best fits the log likelihood function with a flat prior. to the ambiguity curve at different scales of and . These are the scales of interest for signals with strength or (see, Eq. (19)). [Note that Fig. 2 is computed with the reference frequency at 40 Hz. For our results, except for this figure, we choose the reference frequency at 100 Hz. The structure illustrated here is present but less significant at 100 Hz; see appendix B.]

The shape of the ambiguity function has structure on multiple scales. The neighborhood of suggests a much sharper peak than what is seen at the scale. A Fisher matrix computed from formal waveform parameter derivatives (defined in the limit ) or finite difference such that can be overly optimistic about how well can be measured for a signal with .

Therefore, we wish to define an effective Fisher matrix from the curvature of the ambiguity function on the scales of interest. For example, in the case of two parameters, the fitting function is

 P∗=Pmax+p1δλ21+p2δλ22+p12δλ1δλ2 , (25)

where and are fitting coefficients and . We calculate the effective Fisher matrix as

 (^Γij)eff=−∂2P∗∂λi∂λj . (26)

In some cases, especially when a parameter is poorly determined, the variation of the ambiguity function with a parameter may not be well-described by a quadratic. See, for example, the top panel of Fig. 5. Therefore, we find it useful to employ an “iterative” procedure to find the parameters that are amenable to a quadratic fit. For each parameter , we compute the one-dimensional curve and fit it against . This determines the diagonal elements of the effective Fisher matrix, . We discard any parameters that are poorly fit by this quadratic (checked either “by eye” or a with quantitative threshold on the goodness of fit). For the well-fit parameters, we determine the off-diagonal elements of the Fisher matrix by computing the two-dimensional surface and fitting it to a function of the form like Eq. (25) while using the values from the one-dimensional fits, i.e. . We note this method is used primarily as a sanity check to identify any parameters which induce obviously non-quadratic variations in .

Once we have identified the space of all reasonably quadratic parameters, the ambiguity function on that space can be approximated as

 P∗(λ0,λ0+δλ)=1−(^Γij)eff δλi δλj/2 . (27)

Rather than using the iterative procedure described above to find the elements of one at a time, it is straightforward to use a standard least-squares fitting technique to simultaneously solve for all of the . This will also give a better global approximation to the ambiguity function than the iterative approach. As an example of the small yet noticeable differences between these two procedures, Table 6 compares the results for computing the effective Fisher matrix from both “iterative” and “simultaneous” fits to the ambiguity function. Everywhere else in this work (Tables 3, 4, 5, 6, 7, and 9) the effective Fisher matrix is computed by simultaneously fitting all coefficients.

In the cases where is not well-described by a quadratic (e.g., see Fig. 5), we can adopt more complicated expressions to characterize the functional dependence of when these parameters are varied, both in isolation and in correlation with well-constrained variables. As a concrete example, in the absence of higher harmonics the line of sight from the binary is both weakly constrained by observations and nearly separable in from other degrees of freedom333Approximate separability of the line of sight from other degrees of freedom follows only in our well-chosen coordinates, where the binary configuration is specified at . . Specifically, ignoring maximization in time and phase, the overlap between any two lines of sight can be well-approximated by Eq. (B1c) from Osh12 ():

 Pangles ≃ |e2iϕY2(θ)Y2(θ′)+e−2iϕY−2(θ)Y−2(θ′)|√(|Y2(θ)|2+|Y−2(θ)|2)(|Y2(θ′)|2+|Y−2(θ′)|2)

where we use the shorthand and similarly for to reduce superfluous subscripts and where we factor out the common from . This function has wide, nearly flat extrema in for each fixed . On the other hand, in the absence of higher harmonics the line of sight has little impact on the waveform phase versus time away from the orbital plane. We can therefore approximate the ambiguity function for in the top panel of Fig. 5 by

 P≃Pangles[1−12(^Γab)effδλaδλb]−ΓaNδλaδλN , (29)

where the index varies over the line-of-sight parameters and the indices vary over the other parameters and for . With higher harmonics, the functional form above [Eq. (29)] is weakly perturbed by additional angular terms of the form

 P(λ0,λ) ≃Pangles[1−12(^Γab)effδλaδλb]−ΓaNδλaδλN −12(1−cosι)Gϕϕ(ϕ−ϕ0)2 (30) −12Gcc(cosι−cosι20 −Gϕc(ϕ−ϕ0)(cosι−cosι0)

where is a matrix with . This approximation both factors out the leading-order angular dependence and adds additional angular terms with parameter-dependent coefficients, designed to correctly reproduce a -independent result when . Although these terms allow us to correctly reproduce the non-ellipsoidal contours seen in the bottom panel of Fig. 5, Tables 6 and 7 show that this complicated structure only marginally improves the overall fit compared to a purely quadratic approximation . Fit parameters for this more complicated functional dependence are not presented here.

We also considered a more general fit, treating as a parameter. While this parameterization has a significant aesthetic advantage – its effective Fisher matrix is roughly scale-independent when and agrees with analytic calculations – it systematically underestimates in the neighborhood of the maximum. As both fits work well globally, we favor the simpler procedure and adopt except for Table 5.

iii.5 Comparing to standard Fisher matrix results

Despite subtle differences associated with time domain versus frequency domain waveforms, the complex overlap, higher harmonics, and the line of sight, our results for the effective Fisher matrix are directly comparable to earlier results calculated with the stationary phase approximation Poi95 (). For example, for emission along the axis, both the real and complex strain have the form444The two differ for : the complex strain has , while the real strain has . for . As a result, the Fisher matrix can be well-approximated by an identical average over frequency:

 ^Γij = ∫∞−∞df|~h|2(∂aΨ)(∂bΨ)/Sh∫∞−∞df|~h|2/Sh (31)

where we neglect derivatives as small compared to the leading-order phase dependence. Thus, each component of our Fisher matrix must resemble previous results. In fact, as the general definition [Eq. (21)] suggests, each component of (unlike ) depends only on the local response to changing two parameters , no matter how many parameters exist in the model. Therefore, the Fisher matrix for an identical model with more parameters will have, as a submatrix, the Fisher matrix for the smaller model. By contrast, other methods for expressing uncertainty like the covariance depend simultaneously on all terms in . For low-mass binaries, the Fisher matrix is well-known to be poorly conditioned, with eigenvalues spanning several orders of magnitude. We therefore preferentially compare the component-by-component Fisher matrix, rather than the covariances , when comparing results. When presenting results, we provide several significant figures to insure all eigenvalues of remain positive-definite.

Our complex overlap maximizes over time and polarization. The analytic Fisher matrix calculated from the stationary phase approximation [Eq. (31)] has time and phase as parameters. To account for maximizing over those parameters, we transform the full Fisher matrix to a smaller-dimensional matrix which projects out those dimensions:

 (^Γab)max = ^Γab−^ΓaCQCD^ΓCb (32) QCD ≡ [^ΓCD]−1 (33)

where run over the variables and all other variables. In these expressions, the matrix is the inverse of the projection of into the subspace.

iii.6 Comparing to posteriors

Standard parameter estimation techniques like Markov-Chain Monte Carlo produce samples of the full posterior probability distributions, including postprocessed data products like one-dimensional standard deviations and two-dimensional covariances . For strong signals with well-isolated probability distributions, our one-dimensional standard deviations and covariances are directly comparable, for identical binaries. For fainter signals with broad probability distributions, our results will describe part of the posterior, in the neighborhood of one extremum.

For brevity, we have explicitly eliminated two parameters – event time and polarization – and make no predictions about any correlation including them. We will revisit these parameters, along with asymmetric detector response, in a subsequent publication.

iii.7 Numerical and systematic effects

At the very smallest scales, delicate implementation-dependent choices can also introduce artificial structure into the ambiguity function. We have already extensively described how the choice of reference frequency introduces (coordinate-dependent) structure. Less physically, the sampling rate for the waveforms can produce artificial small-scale structure; to avoid this effect we sample at a variety of data rates, typically either (for ) or (for ). Finally, the ambiguity function can also be impacted by our choice for the starting and ending frequency. For our calculations, we start integrating the waveform and integrate over all power above . In our experience, this procedure best mimics the real data processing used in initial LIGO searches. However, a not-insignificant amount of power is present between and ; if included in the integral, the overlap would differ by , comparable to some fine-scale structures of interest. At the other extreme, we terminate our evolution at the minimum-energy circular orbit (MECO), where the binary energy ceases to decrease monotonically.

One small, subtle, but important effect is the nonzero overlap of the waveforms along the axis555Such waveforms would have a real overlap of unity, but the complex overlap is expected to be zero. See Appendix A.. For example, for our non-precessing binary, we find that for otherwise identical parameters

 maxtc,ψ|⟨h(^z)|h(−^z|tc,ψ)⟩|≃1.7×10−3 (34)

Equivalently, in the language of single-detector real overlaps, the sine and cosine chirps are not precisely orthogonal, for the same orbital phase666The real overlap is usually performed in conjunction with an explicit maximization over time and orbital phase, for a single harmonic. For that situation, this subtlety does not occur.. To a first approximation the signals are basis waveforms; the waveform along any line of sight is a superposition of the two. Because these two signals are not orthogonal, the ambiguity function generally has fine-scale structure with , associated with the overlap of these two directions. For example, on this scale and below, the overlap between two non-precessing waveforms with just emission is no longer well-described by Eq. (III.4). Instead, the ambiguity function gains additional fine-scale structure in angle. While extremely useful, for our purposes this result means that on sufficiently small scales , the complex overlap will have additional structure compared to “conventional” investigations of single-detector, optimally-oriented overlaps (i.e., overlaps of two real signals, extracted along ). In particular, this nonzero overlap is partially responsible for the small-scale structure seen in Fig. 3.

Because of the many subtle interpretation and implementation issues associated with the smallest ambiguity scales, while we investigate the value of effective fitting to fine scales (e.g., ), for simplicity we emphasize results for the scale relevant to most detection events ().

Iv Results

Using a small set of fiducial simulations, we compare non-precessing and precessing signals against their immediate neighbors, mapping out an ambiguity function in each -dimensional parameter space.

Higher harmonics perturb the ambiguity function by a quantifiable amount (i.e., for depending on and the harmonic). As has previously been shown elsewhere, we find that higher harmonics break degeneracies present in non-precessing, leading-order signals Bro07 (); Lan06 (); Por08 ().

As described below, we generally find small but significant scale-dependent disagreement with the conventional stationary-phase Fisher matrix calculation, even in the absence of higher harmonics. Motivated by Figs. 2 and 3, as well as the Appendix and Figs. 10 and 11, we suspect that most scale dependence is introduced by suboptimal coordinates and can be minimized by a better choice of reference frequency. Despite our best attempts to find coordinates well-adapted to the problem, the change in going from to for leading-order waveforms is comparable to the change in going from leading-order waveforms to higher harmonics.

iv.1 Zero spin

For a system without spin, higher harmonics principally provide information about the line of sight. For clarity, we will first discuss the most immediately relevant scale (). Comparing the solid (without higher harmonics) and dotted (with higher harmonics) curves on the top panel of Fig. 4, we immediately see that higher harmonics provide little new information about intrinsic parameters, all other things being equal. Equivalently, looking at Table 3, the effective Fisher matrix on the two-dimensional parameters without and with higher harmonics are similar to each other, as well as to a standard Fisher matrix computed using stationary phase approximation waveforms (labeled as ). By contrast, as illustrated by the dramatic difference between the top and bottom panel in Fig. 5, higher harmonics produce a dramatic qualitative change in how well the line of sight can be measured.

Both with and without higher harmonics [Fig. 5], the line of sight is very difficult to measure, particularly at the expected relevant scale (i.e., SNR of around ). In both cases, the ambiguity function has a broad, extended, asymmetric extremum. In the absence of higher harmonics, the ambiguity function cannot be usefully described by a locally quadratic approximation, even a effective one. Nonetheless, by understanding the expected dependence on angle (and by adopting coordinates in band), we can propose a physically-well-motivated fitting function, both for the purely angular dependence and for the correlations between line of sight and other parameters [Eq. (30)]. This fitting function works extremely well when higher harmonics are neglected. When higher harmonics are included, a quadratic approximation sufficies, as the local extremum is much narrower. In the latter case, Table 4 provides the fitting parameters needed to reconstruct the full multidimensional fit. In fact, our well-chosen reference frequency produces a nearly separable fit, with zero off-diagonal terms [e.g., in Eq. (30)]. We will return to this simple structure frequently below.

Several effects besides higher harmonics also introduce fine-scale structure into the ambiguity function. As demonstrated in Fig. 2, the choice of reference frequency can introduce strong, scale-dependent features into the ambiguity function. We chose a reference frequency at to reduce its effect, but have not eliminated it completely. The nonzero overlap between the modes is another such effect. Hence, we are not surprised that our effective fitting parameters change as we reduce the range of used in the fit, even in the absence of harmonics; see Table 3.

The effect of scale dependence is fairly mild: the eigendirections for and agree, only the eigenvalue scale changes. For comparison, we also considered an alternate fitting technique that allowed the single best fit point to have ; see Fig. 3, and Table 5. While this method leads to aesthetically pleasing results similar to analytic calculations, this fit systematically underestimates near the maximum-likelihood point and does not completely eliminate the trend towards different fitting parameters on the smallest scales. We henceforth adopt .

To facilitate approximate comparisons with prior work, Tables 3 and 4 provide one-dimensional standard deviations and correlation coefficients . Unless otherwise stated, these quantities are derived solely from the Fisher matrix fits from that same table. For example, in Table 3, the “measurement errors” follow from inverting the matrix shown, while in Table 4 they follow from inverting a full matrix.

iv.2 Aligned spin

The aligned-spin results at cannot be easily compared with the corresponding zero-spin results (). On the one hand, the component-by-component Fisher matrix coefficients like will differ significantly, as the waveform phasing changes as a function of and hence so does . On the other hand, the aligned-spin results allow a new parameter (spin) that was treated as fixed for the case, with nontrivial coupling to the other intrinsic parameters. The one-parameter uncertainties are dramatically increased by including this previously-neglected systematic effect.

Both at leading and higher-order, our effective fit to the ambiguity function is complicated by the wide range of scales in , even for fixed line of sight. As is well-known from previous Fisher matrix calculations with aligned-spins Poi95 (); Cut94 (), the ambiguity function in has strong correlations, producing a narrow and extended extremum. For the specific example described by our effective Fisher matrix, the submatrix has eigenvalues , describing a strong hierarchy of scales. For our purposes, Fig. 6 demonstrates our fiducial aligned-spin binary cannot be distinguished from binaries with spin : for each in this range, suitable combinations of exist with high overlap. For these extremely extended ambiguity ellipsoids, a fit that reproduces over the full range in might require a more generic functional form than the one adopted so far: a quadratic with constant coefficients in the neighborhood of the fiducial binary. Effectively speaking, however, these additional degrees of freedom add little information with considerable expense. We will explore more complicated effective dependence in a subsequent publication.

As in the zero-spin case, we find significant differences on the smallest scales in , in a fashion that depends on the reference frequency. Given the number of dimensions, complex functional form, sensitivity to numerical implementation like the sampling rate, and less immediate observational relevance, we defer a detailed discussion of fine-scale effects to a subsequent paper.

iv.3 Precessing spin: Case 1

For the first of our two fiducial precessing binaries, we find higher harmonics provide little added information beyond the constraints produced in the non-precessing case. This unfortunate but expected result can be seen, for example, from the one-dimensional covariances in Table 9; from the effective Fisher matrix coefficients ; or from their comparable sequences of eigenvalues. That said, even in the absence of higher harmonics, the ambiguity function for a precessing binary has simpler structure than the non-precessing result, with reduced correlations among the “intrinsic” parameters; a somewhat less extreme hierarchy of scales (i.e., eigenvalues)888We can always rescale our eigenvalues by rescaling our coordinate units. However, in these units all parameters have a prior range of order unity, and physical meaning. As a result, our eigenvalues also have meaning: the coordinate combinations corresponding to the smallest eigenvalues have minimal impact on the overlap and can be ignored.; and roughly speaking a more quadratic ambiguity function.

In fact, for a precessing binary the previous clear separation between “intrinsic” and “geometric” parameters breaks down. As each instant the opening angle of the precession cone of around depends on the relative magnitude of and , as well as on their (nearly conserved) misalignment angle . The magnitudes of and are essentially intrinsic parameters, characterizing the binary masses and BH spin; therefore, we expect the precession cone opening angle to be intimately correlated with the intrinsic parameters. At the same time, the precession cone opening angle must be intimately connected to the “geometric” parameters that define the orientation of the binary at the reference frequency: . Specifically, the orientation of the precession of relative to the line of sight is characterized by the two angles (setting the orientation of )999The other angle needed to specify the orientation of is equivalent to a rotation around the line of sight, i.e. the polarization angle. As the complex overlap maximizes over this angle, it is explicitly removed as a parameter. and (fixing the orientation of along the precession cone at the reference frequency). The orbital phase at the reference frequency fixes the binary’s geometry in band. This cross-coupling between intrinsic and extrinsic parameters has quantitative consequences for the Fisher matrix. Roughly speaking, the two new eigenvalues introduced into by allowing precession (here with values , associated with the and parameters) can be expected to lie between the very large () and very small () eigenvalues associated with the manifestly intrinsic () and extrinsic () parameters.

For non-precessing binaries, the choice of a Hz reference frequency nearly separated intrinsic and extrinsic parameters. A suitable choice of reference frequency may yet further reduce the off-diagonal terms in our effective Fisher matrix. For the present coordinates, however, we cannot cleanly decompose parameters into “intrinsic” and “geometric” parameters. Table 9 shows correlation coefficients calculated by omitting the (nearly unmeasurable) coordiante in ; no obvious block-diagonal form occurs.

iv.4 Precessing spin: Case 2

By contrast to the relatively simple ambiguity functions seen so far, our second set of binary parameters produces a significantly more complicated ambiguity function, particularly in geometric parameters. For example, Fig. 7 shows the ambiguity function versus and all other parameters fixed, for “case 2”. As in case 1, we have a highly symmetric binary starting with , and the line of sight in the same plane at our reference frequency. However, in this case we start with parallel to the line of sight, rather than perpendicular to it. The ambiguity function shows extreme sensitivity to the initial conditions and highly nonquadratic behavior. These differences occur despite the considerable similarity between case 1 and case 2: the two are, to an excellent approximation, the same configurations, just slightly offset in time. By contrast, the change in ambiguity versus is well-described by a quadratic form.

This extreme scenario demonstrates that even an effective Fisher matrix has limits: sometimes, a more generic functional form including higher-order correlations must be used on relevant scales. That said, highly nonquadratic behavior only occurred for a high-symmetry binary. We expect typical binary initial conditions will produce nearly quadratic ambiguity functions.

V Conclusions

In this paper we have used case studies of a coherent, two-detector ambiguity function to estimate how much and what kind of information higher harmonics provide about BH-NS binaries. Given the high dimension of and severe degeneracies that plague the problem, we perform a tractable, idealized calculation instead of the straightforward but less easily understood explicit source in a real multi-detector network. Specifically, we place a single source directly overhead an idealized detector pair, equally sensitive to both polarizations. For all binaries, we find that higher harmonics provide little additional information about the binary’s intrinsic parameters. Instead, at best they provide information about the orientation of the source relative to the line of sight (e.g., for non-precessing binaries). Notably, higher harmonics make it easier to exclude a non-precessing binary with .

When possible, we estimate the two-point function