Measuring the primordial gravitational-wave background in the presence of astrophysical foregrounds
MMeasuring the primordial gravitational-wave backgroundin the presence of astrophysical foregrounds
Sylvia Biscoveanu,
1, 2, 3, a
Colm Talbot,
4, 2, 3
Eric Thrane,
2, 3 and Rory Smith
2, 3 LIGO, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA School of Physics and Astronomy, Monash University, VIC 3800, Australia OzGrav: The ARC Centre of Excellence for Gravitational-Wave Discovery, Clayton, VIC 3800, Australia LIGO, California Institute of Technology, Pasadena, CA 91125, USA
Primordial gravitational waves are expected to create a stochastic background encoding informa-tion about the early Universe that may not be accessible by other means. However, the primordialbackground is obscured by an astrophysical foreground, consisting of gravitational waves from com-pact binaries. We demonstrate a Bayesian method for estimating the primordial background in thepresence of an astrophysical foreground. Since the background and foreground signal parametersare estimated simultaneously, there is no subtraction step, and therefore we avoid astrophysicalcontamination of the primordial measurement, sometimes referred to as “residuals.” Additionally,since we include the non-Gaussianity of the astrophysical foreground in our model, this methodrepresents the statistically optimal approach to the simultaneous detection of a multi-componentstochastic background.
Introduction. —Detection of a cosmologicalgravitational-wave background from the early Universeis one of the most ambitious goals of gravitational-waveastronomy. Broadly speaking, there are two scenarioswhich may give rise to primordial backgrounds: in-flationary scenarios and phase-transition scenarios [1].Inflationary models in general produce a gravitational-wave background through the amplification of vacuumfluctuations [2–5]. In the simplest form of inflation, thedimensionless energy density of the background,Ω gw ( f ) ≡ ρ c dρ gw d ln f , (1)is expected to be Ω gw ( f ) ≈ − across many ordersof magnitude in frequency f [1]. Here, dρ gw is thegravitational-wave energy density between f and f + df while ρ c is the critical energy density for a flat universe.Such a low value of Ω gw is unlikely to be directlydetected by all but the most ambitious space-basedgravitational-wave detectors [6, 7]. However, in mod-els with either non-standard inflation or non-standardcosmology, it is possible to generate inflationary back-grounds accessible by current detectors [8]. Alterna-tively, it may be possible for the inflaton to decay non-perturbatively through parametric resonance. This pro-cess, known as preheating, may produce a potentiallydetectable gravitational-wave background through ex-plosive particle production, peaking as high as Ω gw ≈ − [9, 10]. In reality, the physics of inflation is highlyuncertain. Indirect detection, via the observation of B -modes in the cosmic microwave background, provides analternative means of observing inflationary gravitationalwaves [11].Phase transitions in the early Universe may producegravitational waves if they are strongly first-order [12–15]. The peak frequency of the gravitational-wave energydensity spectrum f is related to the energy scale of the transition T ∗ [1, 13, 16]: f ≈
170 Hz (cid:18) T ∗ GeV (cid:19) . (2)Thus, the detection of a primordial background from aphase transition by either an audio-band or millihertzgravitational-wave detector, would probe physics at en-ergy scales inaccessible by colliders, corresponding to atime when the Universe was only (cid:38) − s old. The en-ergy density created from phase transitions depends onmodel-dependent details, but numerical simulations andscaling arguments suggest that Ω gw ( f ) ≈ − ± fora strongly first-order transition [17]. This is just belowthe projected sensitivity of advanced detectors operatingat design sensitivity [18, 19], but well within the rangeof space-based detectors and proposed third-generationterrestrial detectors [20–22].Astrophysical foregrounds are interesting in their ownright since they contain valuable information about thepopulation properties of compact binaries at high red-shifts [23–25]. However, they complicate searches for theprimordial background. Recent observations of mergingcompact binaries [26–31] by the Advanced LIGO [18]and Virgo [19] detectors imply that primordial back-grounds are masked by much brighter astrophysical fore-grounds [32, 33]. Binary black holes and binary neutronstars each produce astrophysical foregrounds of Ω gw ( f =25 Hz) ≈ − with α = 2 / resolvable with currentdetectors, meaning that some of the events contributingto the foreground are unambiguously detectable. As de-tectors improve, a greater fraction of the foreground isresolved. The most ambitious proposed detectors willresolve essentially every compact binary in the visibleUniverse [6, 7, 38]. These astrophysical foregrounds arenon-Gaussian because the signals do not combine to cre-ate a random signal, characterized only by its statisti- a r X i v : . [ a s t r o - ph . H E ] S e p cal properties. Rather, binary black holes merge every2 −
10 minutes while binary neutron stars merge every4 −
62 s [33]. While there are likely to be many binaryneutron star signals in the LIGO/Virgo band at any giventime, they are nonetheless distinguishable based on theirdifferent coalescence times [38, 39].Previous proposals to disentangle the primordial back-ground from the astrophysical foreground utilize the con-cept of subtraction. The idea, pioneered in [40], is tomeasure the parameters of each resolved compact binaryin order to subtract the gravitational waveform from thedata. Inevitably, this results in “residuals,” or systematicerror from imperfect subtraction. However, the residualscan be “projected out” using a Fisher matrix formalism,which reduces the level of contamination [41]. While [40]considered the subtraction problem in the context of theambitious Big Bang Observer, more recent work has ex-plored the possibility of carrying out subtraction usingthird-generation detectors that may be built in the not-so-distant future [38, 42].One limitation of the subtraction paradigm is thatone can only subtract binary signals that are clearly re-solved. Weak, unresolved signals are not subtracted andtherefore contaminate the measurement of the primor-dial background, introducing a systematic error. Whilebinary black hole mergers will be more easily resolv-able, sub-threshold binary neutron star mergers willimpact the sensitivity of third-generation ground-basedgravitational-wave detectors to the primordial back-ground [43].Other analyses have proposed methods for the simul-taneous measurement of stochastic backgrounds with dif-ferent spectral shapes [44, 45]. However, none of thesemethods account for the non-Gaussianity of the astro-physical foreground, resulting in a decrease in the sensi-tivity of the search. Here, we present a Bayesian formulation in which theprimordial background and the astrophysical foregroundare measured simultaneously. Our method estimates theastrophysical foreground from both resolved and unre-solved binaries, which ensures that our measurement ofthe primordial background is free from bias. The methodcan therefore also include the contributions from highsignal-to-noise ratio compact binaries. Because we modelthe non-Gaussianity of the astrophysical foreground as in[46], this method provides a unified, statistically optimalapproach to the simultaneous detection of compact bina-ries and the primordial background.
Formalism. —We seek to measure a cosmologicalstochastic background described by two power-law pa-rameters: Ω gw ( f ) = Ω α (cid:18) f
25 Hz (cid:19) α . (3)Here, α is a power-law index while Ω α is the amplitude.The background is obscured by a foreground of mergingcompact binaries. Each binary is described by a vector offifteen parameters θ including properties such as the com-ponent masses and the sky location. We only considerbinary black hole mergers for this analysis and assumethat the population distribution for the binary param-eters π ( θ ) is known, but later discuss how the methodcan be generalized to relax these assumptions. Since wewant our formalism to include sub-threshold events, thenumber of compact binary signals in the data is, by as-sumption, unknown.Following [47] and [46], the likelihood of observingfrequency-domain strain data, s i,k , with a Gaussianstochastic background characterized by the parameters(Ω α , α ) and a compact binary coalescence with signal h k ( θ ) is derived by marginalizing over the random Gaus-sian strain perturbation of the background: L ( s i,k | θ, Ω α , α ) = 1det( πT C k (Ω α , α ) /
2) exp (cid:18) − T ( s i,k − h k ( θ )) † C − k (Ω α , α ) ( s i,k − h k ( θ )) (cid:19) , (4)Here, we assume that the data is divided into segmentsof duration T labeled with index i . The frequency depen-dence is denoted with the index k such that s i,k = s i ( f k ).The strain data in each segment, s i,k , and the binarysignal model, h k ( θ ), are vectors with one entry for eachdetector in some network: s i,k = (cid:32) s (1) i,k s (2) i,k (cid:33) h k ( θ ) = (cid:32) h (1) k ( θ ) h (2) k ( θ ) (cid:33) . (5)We use a network consisting of two interferometers inthis analysis, but the framework presented here can be extended to include more detectors. The frequency-dependent covariance matrix, C k , includes contributionsfrom both the detector noise power spectral density(PSD) P I ( f k ) and the the primordial background energydensity: C k = (cid:18) P ( f k ) + κ ( f k )Ω gw κ ( f k )Ω gw κ ( f k )Ω gw P ( f k ) + κ ( f k )Ω gw (cid:19) , (6)where κ IJ ( f k ) ≡ γ IJ ( f k ) 3 H π f (7)converts the primordial background energy density Ω gw into a (signal) strain power spectral density [48, 49]. Thevariable γ IJ ( f k ) is the overlap reduction function for de-tector pair IJ , encoding the geometry of the detectornetwork [47, 50, 51]. It is normalized to γ II = 1 for coin-cident and coaligned detectors with perpendicular arms.Additionally, H is the Hubble constant. Combining datafrom many frequency bins, the likelihood is the productof the individual-frequency likelihoods: L ( s i | θ, Ω α , α ) = m (cid:89) k L ( s i,k | θ, Ω α , α ) . (8)For an astrophysical non-Gaussian foreground, we areinterested in determining the fraction of segments con-taining a signal, ξ , rather than the binary parameters, θ , for a particular segment. We call ξ the “duty cy-cle” following [46]. To be precise, we say that segment“contains” a binary signal if the time of coalescence fallsinside the segment. (We discuss later how this definitionextends to cases where there are multiple binary signalsin the frequency band of the instrument simultaneously.)In this case, the likelihood in Eq. 4 can be marginalizedover the binary parameters θ to obtain L ( s i | Ω α , α, ξ ) = ξ L S ( s i | Ω α , α ) + (1 − ξ ) L N ( s i | Ω α , α ) , (9)where we have defined the marginalized signal and“noise” likelihoods as: L S ( s i | Ω α , α ) = (cid:90) dθ L ( s i | θ, Ω α , α ) π ( θ ) (10) L N ( s i | Ω α , α ) = L ( s i | θ = 0 , Ω α , α ) . (11)The θ = 0 appearing in the expression for L N indi-cates that the noise likelihood is functionally identicalto the signal likelihood if we set the compact binarysignal, h k ( θ ), equal to zero. We note that the name“noise likelihood” is somewhat of a misnomer: while thelikelihood assumes no binary signal, it does allow fora Gaussian background signal. Readers should under-stand the phrase “noise likelihood” to refer to noise +a low-level Gaussian stochastic background. We assumethat the probability of observing one binary black holemerger event in a single segment is much less than one,so that the probability of observing two events is negli-gibly small, which is a reasonable assumption for BBHmergers [33, 38, 46]. We discuss how this assumption canbe relaxed later.For an ensemble of N data segments, { s } , the totallikelihood is given by multiplying the likelihoods for in- dividual segments: L ( { s }| Ω α , α, ξ ) = N (cid:89) i L ( s i | Ω α , α, ξ ) . (12)This joint likelihood function for (Ω α , α, ξ ) defined inEq. 12 is the product of N single-segment likelihoodfunctions, each of which contains a signal sub-hypothesis(with probability ξ ) and a noise sub-hypothesis (withprobability 1 − ξ ).To obtain joint posteriors on (Ω α , α, ξ ), we apply Bayestheorem: p (Ω α , α, ξ |{ s } ) = π (Ω α , α, ξ ) Z L ( { s }| Ω α , α, ξ ) , (13)where Z is the Bayesian evidence given by marginalizingthe total likelihood over the stochastic parameters, Z = (cid:90) d Ω α dα dξ L ( { s }| Ω α , α, ξ ) π (Ω α , α, ξ ) , (14)and π (Ω α , α, ξ ) is the prior. Demonstration. — We demonstrate this formalism withmock data. Assuming a two-detector network of theLIGO Hanford and Livingston observatories operatingat design sensitivity [18], we simulate data for 101 seg-ments each with a duration of 4 s. Each segment con-tains uncorrelated Gaussian noise colored by the noisePSD P ( f k ) of the interferometers as well as correlatedGaussian noise colored by the signal power spectral den-sity of the primordial background. The correlated noiseis simulated such that the cross power spectral density isgiven by κ IJ ( f k )Ω gw ( f k ) for a cosmological backgroundcharacterized by (log Ω α = − , α = 0), where we uselog ≡ log throughout. While this amplitude is severalorders of magnitude higher than that expected for pri-mordial backgrounds, we have chosen this value so thatour simulated cosmological signal corresponds to an un-ambiguous primordial-background detection with signal-to-noise ratio (SNR) of ∼ . α = 0 istypical for cosmological backgrounds [52].Next, we randomly assign binary black hole mergersto 11 of our simulated segments for a corresponding dutycycle of ξ = 11 /
101 = 0 .
11. This duty cycle is higherthan would be expected based on the current estimatesof the BBH merger rate [33], but is just chosen for thepurposes of our demonstration. The chirp mass, M = ( m m ) / ( m + m ) / , (15)is drawn from a uniform prior over the range (13 , M (cid:12) .The prior for the symmetric mass ratio, η = m m ( m + m ) , (16)is uniform over (0 . , . ∝ d L between 500and 5000 Mpc. This results in a range of network opti-mal SNRs between 2.06 and 12.17 with a median of 3.54.Only the signal with the highest SNR corresponds to aconfident detection (with network optimal SNR > < P ( f k )+ κ II ( f k )Ω gw . This results in a decrease in the sen-sitivity of the search, which we mimic in our demonstra-tion by fixing the diagonal terms to the sum of the knownnoise PSD and the signal power from the simulated cos-mological background. Hence, the diagonal terms do notcontribute to the estimation of the (Ω α , α ) parameters.Evaluating the likelihood in Eq. 12 poses a computa-tional challenge due to the product over N single-segmentlikelihoods. To overcome this issue, we use likelihoodreweighting [53] to evaluate the marginalized signal like-lihood (Eq. 10) and the noise likelihood (Eq. 11) on agrid in (Ω α , α ). For each segment we use the cpnest [54]nested sampler as implemented in the Bilby [55, 56]package to obtain posterior samples for the binary pa-rameters using the likelihood in Eq. 4 under the as-sumption that there is no Gaussian background present:Ω α = 0. The priors for the binary parameters are thesame as those used to generate the BBH injections previ-ously described. We use the IMRPhenomPv2 waveformmodel [57–59] for the compact binary signal, h k ( θ ), inboth the simulation and recovery.The marginalized signal likelihood for each segment ata particular value of (Ω α , α ) is calculated via a MonteCarlo integral over the n posterior samples obtained inthe original sampling step: L S ( s i | Ω α , α ) = Z ,i n n (cid:88) j L ( s i | θ j , Ω α , α ) L ( s i | θ j , Ω α = 0) , (17)where Z ,i is the evidence calculated by the sampler usingthe likelihood where Ω α = 0: Z ,i = (cid:90) dθ L ( s i | θ, Ω α = 0) π ( θ ) (18)In the language of importance sampling, the Ω α = 0likelihood is the “proposal” while the Ω α > α , α, ξ ), with the orange lines showing the true valuesused in the simulated data. The 90% credible region isshown in light blue and the 50% credible region in darkblue.be directly evaluated on the same grid in (Ω α , α ) as thereweighted signal likelihood. We use a 50 ×
50 grid rang-ing from log Ω α ∈ [ − , −
4] and α ∈ [0 , ξ ,with 100 values ranging from [0,1]. The full likelihood inEq. 12 is then calculated by multiplying the individualthree-dimensional grids from each data segment. Fig. 1shows the marginalized likelihoods for the cosmologicalbackground parameters (Ω α , α ) as well as ξ obtained us-ing all 101 simulated data segments. We recover val-ues for all three parameters that are consistent with thetrue values used in the simulation: log Ω α = − . +0 . − . , α = 0 . +1 . − . , and ξ = 0 . +0 . − . , where the uncertaintyis the 90% credible interval calculated using the highestprobability density method.In addition to successfully measuring the parameterscharacterizing both the astrophysical foreground and thecosmological stochastic background simultaneously, wealso calculate a Bayes factor comparing the cosmologicalsignal hypothesis to the no-signal hypothesis. This quan-tifies to what extent the model where Ω α = 0 is statisti-cally disfavored compared to the model where (Ω α , α ) cantake on any of the values on our grid. In the high-SNRlimit, the natural log of the Bayes factor is proportionalto the square of the SNR familiar from frequentist cross-correlation searches, ln BF SN ∼ SNR / α andlog Ω α to be uniform across the ranges covered by thegrid. The “noise” evidence is evaluated by integratingEq. 12 assuming that Ω α = 0: Z N = (cid:90) dξ (cid:89) i ξ Z ,i + (1 − ξ ) Z N,i , (19)where Z N,i is the likelihood in Eq. 4 evaluated with both Ω α = 0 and h k ( θ ) = 0. We obtain ln BF SN = ln Z S − ln Z N = 11 .
16, which is consistent with the naive scalingbased on SNR for a signal with SNR = 5 . Discussion. —In this paper we have demonstrateda new method for simultaneously detecting two dis-tinct stochastic gravitational-wave backgrounds—a non-Gaussian astrophysical foreground from sub-thresholdmerging binary black holes and a Gaussian cosmologi-cal background. Our method models both contributionssimultaneously, so that subtraction of the foreground isnot required. Additionally, this is the statistically opti-mal method for detecting a stochastic background con-sisting of both a Gaussian and non-Gaussian component,resulting in significant improvements in the estimatedtime-to-detection of the astrophysical foreground com-pared to other methods for multi-component analyses,as described in [46]. However, in the absence of a non-Gaussian foreground, we find that there is no statisticaladvantage to using the fully Bayesian method comparedto the standard cross-correlation method. Based on thecomparison of the signal-to-noise Bayes factor and SNRfor the presence of the Gaussian background, the twomethods yield a similar level of statistical confidence, tothe extent that it is possible to compare frequentist andBayesian detection statistics.In our demonstration, we assume that the samplingpriors chosen for the binary black hole parameters, π ( θ ),match the true population distribution. In order to avoidbiases that would be introduced due to a mismatch be-tween the population distribution and the sampling prior,our method could be amended to instead measure thesepopulation priors simultaneously with the cosmologicalbackground parameters and duty cycle, following the for-malism described in [61]. This would amount to addingadditional hyper-parameters to the marginalized signallikelihood in Eq. 17: L S ( s i | Λ , Ω α , α ) = Z n n (cid:88) j L ( s i | θ j , Ω α , α ) π ( θ j | Λ) L ( s i | θ j , Ω α = 0 , α ) π ( θ j ) . (20)The hyper-parameters Λ describe the shape of the distri-bution π ( θ | Λ), while the original prior used in the firststep of sampling, π ( θ ), must also be divided out. Thehyper-parameters do not enter the noise likelihood in Eq. 11 because the noise model assumes that each seg-ment contains only the cosmological background with nobinary signal.Evaluating the marginalized signal likelihood in Eq. 20using the same grid-based reweighting technique becomescomputationally prohibitive, since the hyper-parametersΛ drastically increase the dimensionality of the grid. Onepossible solution that has been applied to similar prob-lems in gravitational-wave astronomy could be to builda high-dimensional interpolant [62, 63]. Another promis-ing approach could be to factorize the problem into twoseparate calculations. First, we would carry out popula-tion studies ignoring the stochastic background. SettingΩ = 0 is unlikely to introduce significant bias in our es-timate of the shape of the binary black hole mass, spin,and redshift distributions. Then, once we have the pop-ulation hyper-parameters, we could use the inferred pos-terior predictive distributions for π ( θ | Λ) as priors for theΩ > ∼
15 unresolved binary neutron starsignals in the LIGO band at any given time [33].Because we need to model multiple populations ofmerging binaries simultaneously to avoid contaminationfrom residual power, one possible solution would be totreat the number of binary mergers in a given segmentas a free parameter using a trans-dimensional MarkovChain Monte Carlo algorithm, fitting the binary param-eters for multiple mergers along with the cosmologicalbackground parameters all at once [64]. Another pos-sible method is to work with shorter 0 . ξ presented above. Pre-liminary tests suggests that the presence of other binarysignals, merging at times outside of the segment, havea negligible effect on inferences about the binary merg-ing during the segment. By marginalizing over the bi-nary neutron star parameters in many short segments, itshould be possible to calculate the likelihood of a muchlonger span of data given the stochastic parameters. Weleave investigation of these methods to future work.The formalism we describe and demonstrate as-sumes that the uncorrelated detector noise is Gaussian,while it is known that interferometric gravitational-wavedata suffers from non-Gaussian noise transients calledglitches [66]. This assumption can be relaxed via theintroduction of additional duty cycle parameters to thelikelihood in Eq. 4, characterizing the fraction of seg-ments that contain a glitch in each detector, as describedin [46]. This would increase the computational cost foreach data segment analyzed, but the method is embar-rassingly parallelizable, so the overall wall time for run-ning the analysis does not increase significantly.We also note that limitations in the accuracy of thewaveform model describing the compact binary signalcan leave behind coherent residual power that could biasthe inference of the Gaussian background parameters.Based on current estimates of the uncertainty in numer-ical relativity waveforms [67], this level of contaminationwould likely not affect cosmological backgrounds probe-able with proposed third-generation detectors, but im-provements to waveform modeling would be necessaryto recover unbiased parameter estimates for the weakestbackground models.While we demonstrate our method for the simulta-neous detection of a stochastic background with bothGaussian and non-Gaussian components in the contextof a cosmological background and an astrophysical fore-ground of binary black hole mergers, this formalism canbe applied to any analogous problems. For example, thismethod could be applied to simultaneously measure bothindividual compact binary mergers or a foreground ofthese sources in the frequency band of the space-basedLISA detector [68] on top of the white dwarf confusionnoise background [69, 70]. Our model can also be ex-tended to include multiple Gaussian backgrounds withdifferent spectral shapes through the addition of extraterms in the covariance matrix defined in Eq. 6. One suchexample is the contamination from correlated magneticnoise in a ground-based detector network [50, 71, 72],which has a unique overlap reduction function [73]. Acknowledgements. —The authors thank Jan Harmsand Salvatore Vitale for thoughtful comments on themanuscript. SB and CT acknowledge support of theNational Science Foundation and the LIGO Laboratory.LIGO was constructed by the California Institute ofTechnology and Massachusetts Institute of Technologywith funding from the National Science Foundation andoperates under cooperative agreement PHY-1764464. SBis also supported by the Paul and Daisy Soros Fellowshipfor New Americans, the Australian-American FulbrightCommission, and the NSF Graduate Research Fellow-ship under Grant No. DGE-1122374. ET and RS aresupported by the Australian Research Council (ARC)CE170100004. ET is supported by ARC FT150100281.The authors are grateful for computation resources pro-vided by the OzStar cluster. This article carries LIGODocument Number LIGO-P2000297. a [email protected] [1] M. Maggiore, Phys. Rept. , 283 (2000), arXiv:gr-qc/9909001.[2] A. H. Guth, Phys. Rev. D , 347 (1981).[3] M. S. Turner, Phys. Rev. D , 435 (1997), arXiv:astro-ph/9607066.[4] R. Easther, J. Giblin, John T., and E. A. Lim, Phys.Rev. Lett. , 221301 (2007), arXiv:astro-ph/0612294.[5] M. Guzzetti, N. Bartolo, M. Liguori, and S. Matarrese,Riv. Nuovo Cim. , 399 (2016), arXiv:1605.01615 [astro-ph.CO].[6] J. Crowder and N. J. Cornish, Phys. Rev. D , 083005(2005), arXiv:gr-qc/0506015.[7] S. Kawamura et al. , Class. Quant. Grav. , 094011(2011).[8] V. Mandic and A. Buonanno, Phys. Rev. D , 063008(2006).[9] S. Khlebnikov and I. Tkachev, Phys. Rev. D , 653(1997).[10] L. R. Price and X. Siemens, Phys. Rev. D , 063541(2008).[11] P. Ade et al. (BICEP2, Keck Array), Phys. Rev. Lett. , 221301 (2018), arXiv:1810.05216 [astro-ph.CO].[12] M. Kamionkowski, A. Kosowsky, and M. S. Turner,Phys. Rev. D , 2837 (1994), arXiv:astro-ph/9310044.[13] T. Kahniashvili, A. Kosowsky, G. Gogoberidze, andY. Maravin, Phys. Rev. D , 043003 (2008),arXiv:0806.0293 [astro-ph].[14] C. Caprini, R. Durrer, T. Konstandin, and G. Servant,Phys. Rev. D , 083519 (2009), arXiv:0901.1661 [astro-ph.CO].[15] C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau,P. Schwaller, G. Servant, and D. J. Weir, JCAP ,001 (2016), arXiv:1512.06239 [astro-ph.CO].[16] B. Allen, in Les Houches School of Physics: AstrophysicalSources of Gravitational Radiation (1996) pp. 373–417.[17] J. T. Giblin and E. Thrane, Phys. Rev. D , 107502(2014), arXiv:1410.4779 [gr-qc].[18] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. ,074001 (2015), arXiv:1411.4547 [gr-qc].[19] F. Acernese et al. (VIRGO), Class. Quant. Grav. ,024001 (2015), arXiv:1408.3978 [gr-qc].[20] E. Thrane and J. D. Romano, Phys. Rev. D , 124032(2013), arXiv:1310.5300 [astro-ph.IM].[21] D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Bar-sotti, G. Billingsley, D. A. Brown, Y. Chen, D. Coyne,R. Eisenstein, et al. , in Bulletin of the American Astro-nomical Society , Vol. 51 (2019) p. 35, arXiv:1907.04833[astro-ph.IM].[22] M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Bel-gacem, D. Bertacca, M. A. Bizouard, M. Branchesi,S. Clesse, S. Foffa, J. Garc´ıa-Bellido, et al. , jcap ,050 (2020), arXiv:1912.02622 [astro-ph.CO].[23] V. Mandic, E. Thrane, S. Giampanis, and T. Regimbau,Phys. Rev. Lett. , 171102 (2012).[24] S. Vitale, W. M. Farr, K. Ng, and C. L. Rodriguez, As-trophys. J. Lett. , L1 (2019), arXiv:1808.00901 [astro-ph.HE].[25] T. Callister, M. Fishbach, D. E. Holz, and W. M. Farr,apjl , L32 (2020), arXiv:2003.12152 [astro-ph.HE].[26] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 061102 (2016), arXiv:1602.03837 [gr-qc].[27] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 161101 (2017), arXiv:1710.05832 [gr-qc]. [28] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X , 031040 (2019), arXiv:1811.12907 [astro-ph.HE].[29] B. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J.Lett. , L3 (2020), arXiv:2001.01761 [astro-ph.HE].[30] R. Abbott et al. (LIGO Scientific, Virgo), (2020),arXiv:2004.08342 [astro-ph.HE].[31] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. , L44 (2020), arXiv:2006.12611 [astro-ph.HE].[32] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 131102 (2016), arXiv:1602.03847 [gr-qc].[33] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 091101 (2018), arXiv:1710.05837 [gr-qc].[34] X.-J. Zhu, E. Howell, T. Regimbau, D. Blair, and Z.-H. Zhu, Astrophys. J. , 86 (2011), arXiv:1104.3565[gr-qc].[35] C. Wu, V. Mandic, and T. Regimbau, Phys. Rev. D ,104024 (2012), arXiv:1112.1898 [gr-qc].[36] P. A. Rosado, Phys. Rev. D , 084004 (2011).[37] X.-J. Zhu, E. J. Howell, D. G. Blair, and Z.-H. Zhu, Mon.Not. Roy. Astron. Soc. , 882 (2013), arXiv:1209.0595[gr-qc].[38] T. Regimbau, M. Evans, N. Christensen, E. Katsavouni-dis, B. Sathyaprakash, and S. Vitale, Phys. Rev. Lett. , 151105 (2017), arXiv:1611.08943 [astro-ph.CO].[39] F. H. Vivanco, R. Smith, E. Thrane, and P. D. Lasky,Phys. Rev. D , 043023 (2019).[40] C. Cutler and J. Harms, Phys. Rev. D , 042001 (2006),arXiv:gr-qc/0511092.[41] J. Harms, C. Mahrdt, M. Otto, and M. Priess, Phys.Rev. D , 123010 (2008), arXiv:0803.0226 [gr-qc].[42] A. Sharma and J. Harms, (2020), arXiv:2006.16116 [gr-qc].[43] S. Sachdev, T. Regimbau, and B. Sathyaprakash,(2020), arXiv:2002.05365 [gr-qc].[44] C. Ungarelli and A. Vecchio, Classical and QuantumGravity , S857 (2004).[45] A. Parida, S. Mitra, and S. Jhingan, Journal of Cosmol-ogy and Astroparticle Physics , 024 (2016).[46] R. Smith and E. Thrane, Physical Review X , 021019(2018), arXiv:1712.00688 [gr-qc].[47] J. D. Romano and N. J. Cornish, Living Reviews in Rel-ativity , 2 (2017), arXiv:1608.06889 [gr-qc].[48] B. Allen and J. D. Romano, Phys. Rev. D , 102001(1999).[49] C. M. Mingarelli, S. R. Taylor, B. Sathyaprakash, andW. M. Farr, (2019), arXiv:1911.09745 [gr-qc].[50] N. Christensen, Phys. Rev. D , 5250 (1992).[51] E. E. Flanagan, Phys. Rev. D , 2389 (1993).[52] N. Christensen, Rept. Prog. Phys. , 016903 (2019),arXiv:1811.08797 [gr-qc].[53] E. Payne, C. Talbot, and E. Thrane, Phys. Rev. D , 123017 (2019), arXiv:1905.05477 [astro-ph.IM].[54] W. D. Pozzo and J. Veitch, (2017),DOI:10.5281/zenodo.835874.[55] G. Ashton et al. , Astrophys. J. Suppl. , 27 (2019),arXiv:1811.02042 [astro-ph.IM].[56] I. Romero-Shaw et al. , (2020), arXiv:2006.00714 [astro-ph.IM].[57] M. Hannam, P. Schmidt, A. Boh, L. Haegel, S. Husa,F. Ohme, G. Pratten, and M. Prrer, Phys. Rev. Lett. , 151101 (2014), arXiv:1308.3271 [gr-qc].[58] S. Husa, S. Khan, M. Hannam, M. Prrer, F. Ohme,X. Jimnez Forteza, and A. Boh, Phys. Rev. D , 044006(2016), arXiv:1508.07250 [gr-qc].[59] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Prrer,X. Jimnez Forteza, and A. Boh, Phys. Rev. D , 044007(2016), arXiv:1508.07253 [gr-qc].[60] E. Thrane and C. Talbot, Publ. Astron. Soc. Austral. ,e010 (2019), arXiv:1809.02293 [astro-ph.IM].[61] R. J. E. Smith, C. Talbot, F. HernandezVivanco,and E. Thrane, Mon. Not. Roy. Astron. Soc. ,3281 (2020), https://academic.oup.com/mnras/article-pdf/496/3/3281/33464232/staa1642.pdf.[62] F. Hernandez Vivanco, R. Smith, E. Thrane, P. D. Lasky,C. Talbot, and V. Raymond, Phys. Rev. D , 103009(2019), arXiv:1909.02698 [gr-qc].[63] D. Wysocki, R. O’Shaughnessy, L. Wade, and J. Lange,(2020), arXiv:2001.01747 [gr-qc].[64] P. J. Green, Biometrika , 711 (1995),https://academic.oup.com/biomet/article-pdf/82/4/711/699533/82-4-711.pdf.[65] F. Hernandez Vivanco, R. Smith, E. Thrane, andP. D. Lasky, Phys. Rev. D , 043023 (2019),arXiv:1903.05778 [gr-qc].[66] B. P. Abbott et al. (LIGO Scientific, Virgo), Class.Quant. Grav. , 134001 (2016), arXiv:1602.03844 [gr-qc].[67] M. Prrer and C.-J. Haster, Phys. Rev. Res. , 023151(2020), arXiv:1912.10055 [gr-qc].[68] P. Amaro-Seoane et al. (LISA), (2017), arXiv:1702.00786[astro-ph.IM].[69] P. L. Bender and D. Hils, Classical and Quantum Gravity , 1439 (1997).[70] M. R. Adams and N. J. Cornish, Phys. Rev. D , 022001(2014).[71] E. Thrane, N. Christensen, and R. Schofield, Phys. Rev.D , 123009 (2013), arXiv:1303.2613 [astro-ph.IM].[72] E. Thrane, N. Christensen, R. M. S. Schofield, and A. Ef-fler, Phys. Rev. D , 023013 (2014), arXiv:1406.2367[astro-ph.IM].[73] Y. Himemoto and A. Taruya, Phys. Rev. D96