Parameter estimation for gravitational-wave bursts with the BayesWave pipeline
Bence Bécsy, Peter Raffai, Neil J. Cornish, Reed Essick, Jonah Kanner, Erik Katsavounidis, Tyson B. Littenberg, Margaret Millhouse, Salvatore Vitale
aa r X i v : . [ a s t r o - ph . H E ] A p r Draft version April 11, 2017
Preprint typeset using L A TEX style AASTeX6 v. 1.0
PARAMETER ESTIMATION FOR GRAVITATIONAL-WAVE BURSTS WITH THE BAYESWAVE PIPELINE
Bence B´ecsy , Peter Raffai , Neil J. Cornish , Reed Essick , Jonah Kanner , Erik Katsavounidis , Tyson B.Littenberg , Margaret Millhouse , and Salvatore Vitale Institute of Physics, E¨otv¨os University, 1117 Budapest, Hungary; [email protected] MTA-ELTE EIRSA “Lend¨ulet” Astrophysics Research Group, 1117 Budapest, Hungary Department of Physics, Montana State University, Bozeman, MT 59717, USA Massachusetts Institute of Technology, 185 Albany St, 02139 Cambridge USA LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA NASA Marshall Space Flight Center, Huntsville AL 35812, USA
ABSTRACTWe provide a comprehensive multi-aspect study on the performance of a pipeline used by the LIGO-Virgo Collaboration for estimating parameters of gravitational-wave bursts. We add simulated signalswith four different morphologies (sine-Gaussians, Gaussians, white-noise bursts, and binary black holesignals) to simulated noise samples representing noise of the two Advanced LIGO detectors duringtheir first observing run. We recover them with the BayesWave (BW) pipeline to study its accuracy insky localization, waveform reconstruction, and estimation of model-independent waveform parameters.BW localizes sources with a level of accuracy comparable for all four morphologies, with the medianseparation of actual and estimated sky locations ranging from 25.1 ◦ to 30.3 ◦ . This is a reasonable ac-curacy in the two-detector case, and is comparable to accuracies of other localization methods studiedpreviously. As BW reconstructs generic transient signals with sine-Gaussian wavelets, it is unsur-prising that BW performs the best in reconstructing sine-Gaussian and Gaussian waveforms. BW’saccuracy in waveform reconstruction increases steeply with network signal-to-noise ratio (SNR net ),reaching a 85% and 95% match between the reconstructed and actual waveform below SNR net ≈ net ≈
50, respectively, for all morphologies. BW’s accuracy in estimating central moments ofwaveforms is only limited by statistical errors in the frequency domain, and is affected by systematicerrors too in the time domain as BW cannot reconstruct low-amplitude parts of signals overwhelmedby noise. The figures of merit we introduce can be used in future characterizations of parameterestimation pipelines.
Keywords: gravitational waves – methods: data analysis INTRODUCTIONThe network of Advanced LIGO (aLIGO)gravitational-wave (GW) detectors (Aasi et al. 2015),consisting of aLIGO-Hanford (H1) and aLIGO-Livingston (L1), finished its first observing run (O1) inJanuary 2016. During O1, this network achieved thefirst direct detections of GWs by detecting GW150914(Abbott et al. 2016d) and GW151226 (Abbott et al.2016b), two signals from coalescences of binary blackholes. Besides binary black holes, other astrophysicalsources of GW transients (e.g. core-collapse supernovae,magnetar flares, and cosmic string cusps) are alsotargeted by aLIGO (Abbott et al. 2016g). Searches forgeneric GW transients aim to detect weakly-modeledGW signals (“bursts”) from such systems, as well asfrom binary black holes, and also from as-yet-unknownsources (see e.g. Abbott et al. 2016e). Detections of GW signals will be used to testand constrain models of astrophysical sources (seee.g. Abbott et al. 2016a; Belczynski et al. 2016). Thisusually requires reconstructing the signal waveform fromthe GW detector output and estimating parameters ofthe waveform (see e.g. Abbott et al. 2016f). For sourceswhere an accurate waveform model exists, such as binaryblack holes in circular orbits, this is done by match-ing the detector output with template waveforms (seee.g. Abbott et al. 2016f). In this case, the estimated pa-rameters are astrophysical, e.g., chirp mass and spins.Parameter estimation (PE) for burst signals, where nomodel templates exist, need a different approach. Insuch cases, basis functions are used to reconstruct thewaveform and to estimate model-independent param-eters of it, such as central time and frequency, signalduration and bandwidth. Besides these intrinsic param-eters of the waveform, estimates can also be given on theextrinsic parameters of the source (e.g. sky location).BayesWave (BW) is a pipeline for detecting and char-acterizing GW bursts, that works within the frame-work of Bayesian statistics and uses sine-Gaussianwavelets as basis functions to reconstruct the signal(Cornish & Littenberg 2015). In O1, BW was used asa follow-up PE tool on triggers provided by the coher-ent Waveburst (cWB) search pipeline (Klimenko et al.2008, 2016), which identifies coincident excess power instrain data of multiple GW detectors. Note however,that cWB can also reconstruct the sky location of a GWsource and the waveform of the GW signal, indepen-dently from BW (Klimenko et al. 2011). This providesan opportunity to compare the performances of BW andcWB in PE using the same set of triggers (for the re-sults of this comparison, see Section 3.1). BW is ef-fective in distinguishing GW signals from non-Gaussiannoise artifacts (“glitches”), which enables the combi-nation of the cWB and BW pipelines to achieve high-confidence detections across a range of waveform mor-phologies (Littenberg et al. 2016; Kanner et al. 2016).The estimates of mass parameters and sky location ob-tained by BW for GW150914 have shown to be con-sistent with template-based PE pipelines (Abbott et al.2016e).In this paper we characterize BW’s performance inPE by injecting a large set of simulated signals intosimulated aLIGO noise, and recovering them and theirparameters with BW. The main purpose of this studyis to determine the accuracy of these reconstructionsthat can be achieved with BW. By knowing the accu-racy, future studies can identify the broadest range ofastrophysical models that can be tested with BW, whilefurther improvements of BW can be guided by these re-sults. Among the estimated parameters, we give spe-cial attention to sky location of the GW source, be-cause of its key role in electromagnetic (EM) follow-up observations of GW events (see e.g. Abbott et al.2016c; Singer et al. 2014; Berry et al. 2015; Vitale et al.2017). Sky localization of GW burst sources can be car-ried out with the cWB and LALInferenceBurst (LIB)pipelines (Lynch et al. 2015; Veitch et al. 2015) too. Anextensive analysis on the sky localization performanceof cWB and LIB was published in Essick et al. (2015).Here we present a similar analysis for BW in order tocharacterize its performance and to allow comparisonswith other burst pipelines studied in Essick et al. (2015).Note however, that as we use a reduced set of trig-gers compared to Essick et al. (2015) (for an explana-tion, see Appendix A), our results in Figures 1-4 shouldnot be compared directly with results in Figures 3-6 ofEssick et al. (2015). Instead, to allow direct compar-isons between BW, cWB, and LIB, we repeat our anal- ysis with cWB and with LIB on the same reduced set oftriggers, and present the results in Figure Sets 1-4 (avail-able in the online journal). Also note that new cWBsky localization results for binary black holes presentedrecently (see Vitale et al. 2017) show that cWB’s per-formance has improved significantly for a three-detectornetwork, while it has not changed significantly for thetwo-detector case we present here.We focus on three aspects of BW’s performance: (i)sky localization, (ii) waveform reconstruction, and (iii)estimation of model-independent waveform parameters.In Section 2, we describe the methods used for creatingsimulated signals and noise samples, and used by BW forcarrying out PE. In Section 3, we present results of ouranalyses regarding all (i)-(iii) aspects. We summarizeour findings and highlight some implications in Section4. METHODSWe used software injections to test the PE perfor-mance of BW, i.e. we created mock samples of aLIGOnoise and added simulated GW signals with four differ-ent morphologies to these samples. We then used thesesamples at trigger times provided by cWB as inputs forBW to test what it recovers from the signals embeddedin the mock detector noise. In this section we discussthe characteristics of noise samples and of simulated sig-nals we used (Section 2.1), as well as methods BW usesfor PE (Section 2.2).2.1.
Noise and injections
In this section we summarize characteristics of injec-tions and noise samples we used in our analyses, whichare the same as the ones used in Essick et al. (2015).For further details on this see Section 2, Appendix Cand Table 4 in Essick et al. (2015).In our analysis we considered a two-detector networkconsisting of H1 and L1. We used stationary, Gaus-sian mock noise samples generated using the expected2015 sensitivity curve of aLIGO, thus they have slightlydifferent characteristics than the actual noise collectedduring the O1 run. Projections show that the two LIGOdetectors will operate in the first two months of the sec-ond observing run (O2) with similar sensitivity curvesthey operated with during O1. Thus, we expect thatour results are representative for this first period of O2as well.Our set of software injections consists of signals withfour different morphologies: sine-Gaussians ( SG ); Gaus-sians ( G ); white-noise bursts ( WNB ); and binary blackhole (
BBH ) mergers. This wide range of signal mor-phologies allows us to test the PE performance of BWwith minimal assumptions on the GW signal. The am-plitude distribution of injected signals was chosen such
Table 1 . Number of injected signals for each morphology at different stages ofthe analysis. For details on why BW identified many SG and WNB signals asglitches or Gaussian noise, and how this has been improved for O2, see AppendixA.
SG G WNB BBH
Triggers produced by cWB 1112 256 769 2488Left out to reduce computational costs 0 0 0 -1988Analyzed by BW 1112 256 769 500Identified as glitches or Gaussian noise by BW -779 0 -355 -1
Used in our analysis 333 256 414 499 as to represent a uniform distribution of GW sourcesin volume. Signal injections were distributed uniformlyover the sky and were regularly spaced in time.The number of signals we analyzed was determinedby multiple factors (see Table 1): (i) the BW ver-sion we used runs only on triggers produced by cWB (Abbott et al. 2016e); (ii) we reduced the number of
BBH triggers in order to reduce computational costs;(iii) we only used signals correctly identified as signalsby BW. For details on why BW identified many SG and WNB signals as glitches or Gaussian noise, and how thishas been improved for O2, see Appendix A. SG waveforms are often used to model generic tran-sients (e.g. Abadie et al. 2012), because they are themost localized signals in time-frequency space wheregeneric burst searches (including cWB) operate (seeChatterji 2005). We define SG waveforms with the fol-lowing two equations: h + ( t ) = cos( α ) h rss s f √ πQ (1 + cos(2 φ ) e − Q ) × cos(2 πf ( t − t ) + φ ) e − ( t − t ) /τ (1a) h × ( t ) = sin( α ) h rss s f √ πQ (1 − cos(2 φ ) e − Q ) × sin(2 πf ( t − t ) + φ ) e − ( t − t ) /τ , (1b)where α ∈ [0 , π/
2] is a parameter that sets therelative weights between polarizations h + and h × , h = R ( h + h × )dt is the square of the root-sum-squared strain amplitude chosen as a free parameter inthe amplitude randomization process, f is the centralfrequency, t is the central time, φ is the phase at time t = t , τ is the width of the signal in the time domain,and Q = √ πτ f is the quality factor encoding the char-acteristic number of cycles within duration of the signal. G signals are the special cases of SG signals when f →
0, and are defined as: h + ( t ) = cos( α ) h rss √ τ (cid:18) π (cid:19) / e − ( t − t ) /τ (2a) h × ( t ) = sin( α ) h rss √ τ (cid:18) π (cid:19) / e − ( t − t ) /τ . (2b) Despite their similarity to SG s, these signals pose dif-ferent challenges, because they have their highest am-plitude at f = 0 Hz in frequency domain, and thusthey have most of their power at low frequencies whereaLIGO is less sensitive. WNB waveforms are intended to model a time-localized excess power uniformly distributed in a givenfrequency band, and satisfy: h + , × ( t ) ∝ e − ( t − t τ Z ∞−∞ e − i πft w ( f )d f, (3)where w ( f ) values are randomly drawn from a Gaussianwhite noise within and chosen to be w ( f ) = 0 outsidethe band f ∈ [ f min , f max ]. We generated the right sideof equation (3) independently for the + and × polariza-tions, and normalized them to get h + and h × with thedesired h rss . Unlike signals with the other three mor-phologies, WNB signals are not elliptically polarized,because the procedure used to produce them generates h + and h × independently.The only astrophysical signals we used were BBH swith spins aligned or anti-aligned with the orbital an-gular momentum. We only considered binaries withrelatively high detector-frame total masses ( M tot ∈ [30 , M ⊙ ), because their signals are more compactin time-frequency space, which makes them good tar-gets for generic burst searches. Three different meth-ods have been used for calculating the waveform in thethree different phases of binary evolution: 3.5PN post-Newtonian expansion, numerical relativity, and analyticquasi-normal modes to calculate the inspiral, merger,and ringdown waveforms, respectively (see Ajith et al.2011 and Hannam et al. 2010 for details). 2.2. The BayesWave pipeline
BW uses a trans-dimensional Reversible JumpMarkov Chain Monte Carlo (RJMCMC) algorithm(Green 1995) to explore the following three compet-ing models of the data, and test them with the inputdata samples from each aLIGO detector: i) Gaussiannoise only; ii) Gaussian noise with glitches; iii) Gaus-sian noise with a GW signal. This approach makesBW effective in distinguishing GW signals from glitches(Littenberg et al. 2016), but it also makes BW compu-tationally expensive, and thus in O1 BW was used tofollow-up candidate events from cWB.BW assumes that all signals are elliptically polarized,i.e. h × = ǫh + e iπ/ , where ǫ ∈ [0 ,
1] is the ellipticity pa-rameter, which is 0 for linearly polarized signals and 1for circularly polarized ones. This is a valid assump-tion for many expected astrophysical signals, but notfor our injections with
WNB morphology (see Section2.1). However, for a LIGO-only network, it is often thecase that only a single combination of the two polariza-tions, rather than the separate + and × components,will be detectable, making the elliptical constraint a fairapproximation for many cases.We used the BW version which had been used for theoffline analysis of O1 data to attain a characterizationof BW’s performance during O1, and to support a faircomparison with the versions of other PE pipelines char-acterized in Essick et al. (2015). PE pipelines used bythe LIGO-Virgo Collaboration (including BW) have un-dergone improvements since the beginning of O1 (someof which were motivated by this study). RESULTSIn this section we show how BW performed in differentaspects of PE. These aspects are sky localization (seeSection 3.1), waveform reconstruction (see Section 3.2),and point estimates of waveform central moments (seeSection 3.3).Even though BW’s current (O2) version is more effi-cient in identifying signals (see Appendix A), we usedthe version of BW used during O1 in order to char-acterize BW’s performance during O1, and to allow acomparison of our results with the ones presented inEssick et al. (2015). We only analyzed signals that wereproperly identified as signals by BW (see Table 1). Wepresent a reproduction of results of Essick et al. (2015)for the subset of events we used in this study, to enablea fair comparison of sky localization results (see FigureSets 1-4.).Results presented here depend on the parameter dis-tributions of injected signals defined in Table 4 ofEssick et al. (2015), and on the corresponding detec-tion efficiencies of the combination of cWB and BW pipelines for the different parameter sets. Results areparticularly dependent on the chosen h rss distributionof injected signals, and thus on the network signal-to-noise ratio (SNR net ) distribution of them (see inset ofFigure 5). However, the h rss distribution we chose forthis study is a good approximation for generic burst sig-nals uniformly distributed in volume (see Appendix Cin Essick et al. 2015).3.1. Sky localization
BW computes a skymap defined as the posterior prob-ability density function of the GW source location ex-pressed as a function of celestial coordinates α (rightascension) and δ (declination), denoted by p sky ( α, δ ).Example skymaps for each morphology are shown in Ap-pendix B. Skymaps for all the injections can be foundin the Burst First2Years sky localization Open Data re-lease . There are many possible quantitative measuresfor the “goodness” of source localization, here we imple-ment the ones defined in Essick et al. (2015), i.e. angu-lar offset , searched area , extent and fragmentation . Wereproduced results of Essick et al. (2015) for LIB andcWB using the same subset of events we used in thisstudy (the ones identified as signals by BW) to enablea direct comparison of the results (see Figure Sets 1-4).The first measure is the angular offset ( δθ ), which isthe angular distance between the maximum of p sky andthe true location of the injected signal. Figure 1 showsnormalized histograms of cos( δθ ) for all injections, withthe upper axis showing the corresponding δθ values. Thedistribution has a peak at cos( δθ ) = 1, which suggeststhat BW tends to reconstruct the most probable loca-tion of the source close to the actual source location.There is also a smaller peak at cos( δθ ) = −
1, which in-dicates that it is more likely that BW reconstructs theopposite direction of the sky compared to the locationof the injected signal than a direction perpendicular tothe injected signal’s location. This is due to the factthat opposite directions cannot be distinguished usingthe network antenna pattern which has the same valueat opposite directions because of the near co-alignmentof H1 and L1 detectors (Singer et al. 2014). However,the peak at cos( δθ ) = − δθ ) = 1 because opposite directions are only allowedby the triangulation ring when the source is right above(or below) the detectors, and thus the triangulation ringis a great circle on the celestial sphere. Note that the dis-tributions for different morphologies are very similar toeach other, which means that the angular offset dependsweakly on signal morphology. We show summary statis-tics of δθ distributions for all morphologies in Table 2. Figure 1 . Normalized histograms of angular offsets ( δθ ) forinjections with four different morphologies ( SG , G , WNB , BBH ). Most of the injected signals have cos( δθ ) = 1, whichindicates that BW tends to place the most probable locationclose to the true location. Note that the distributions fordifferent morphologies are very similar to each other, whichmeans that the angular offset does not depend strongly onsignal morphology. The complete figure set (3 figures) show-ing the same plot for cWB and LIB pipelines is available inthe online journal. It is clearly visible that BW performs best for
BBH sig-nals, while SG , G and WNB signals show slightly larger δθ values. Statistical errors on reported values are in theorder of a few percent. Figure Set 1 shows normalizedhistograms of cos( δθ ) obtained with the cWB and LIBpipelines on the subset of signals identified as signals byBW.EM follow-up observations tend to target the point ofthe sky with the highest p sky value first, and continuewith points having lower p sky values. This motivatesthe introduction of the searched area ( A ) as a secondmeasure, which is the total sky area observed beforeaiming a hypothetical telescope at the true location ofthe source: A = Z H ( p sky ( α, δ ) − p ) dΩ (4)where H is the Heaviside step function, p is the valueof p sky at the true location of the source, and dΩ =cos δ d δ d α .We show the cumulative histogram of A for all injec-tions in Figure 2. Histograms for different morphologiesfollow a similar trend, but the curves are shifted alongthe horizontal axis. This can be quantified e.g. with me-dian searched area, which is 252.8 deg for G , 151.0 deg Figure 2 . Cumulative histograms of searched area ( A ). His-tograms for different morphologies follow a similar trend, ex-cept that the curves are shifted along the horizontal axis. Areference curve labeled with SG (LIB) shows results for theLIB pipeline on the subset of SG signals identified as sig-nals by BW. The complete figure set (3 figures) showing thesame plot for cWB and LIB pipelines is available in the on-line journal. for WNB , 121.3 deg for SG , and 99.2 deg for BBH sig-nals. Another difference between morphologies is thatthere is a fraction of
WNB signals with searched areaequal to the whole sky (
A ≃ · deg ). This is due tothe fact that p = 0 for these signals, i.e. the posteriordistribution has no support at the true location of thesource. There are no such signals with SG , G and BBH morphologies. A reference curve labeled with SG (LIB)shows results for the LIB pipeline on the subset of SG signals identified as signals by BW. Note that LIB uses asingle sine-Gaussian to reconstruct the signal, so for SG injections LIB becomes a matched-filtering analysis forwhich better performance is expected, while BW some-times uses more than one sine-Gaussian, because it fa-vors more complex signals. It shows that LIB performedsimilarly, but slightly better for SG signals. We showsummary statistics of A distributions for all morpholo-gies in Table 2. It is clearly visible that BW performsbest for BBH signals, while SG , G and WNB signalsshow significantly larger A values. Statistical errors onreported values are in the order of a few percent. FigureSet 2 shows normalized histograms of A obtained withthe cWB and LIB pipelines on the subset of signals iden-tified as signals by BW.Even if δθ and A are small, the favored sky positionscan still be either well localized or spread out over vari- Figure 3 . Normailzed histograms of the extent ( δθ inj ) ofskymaps for the four different injection morphologies. Theshown distributions are bimodal for all morphologies withpeaks at cos( δθ inj ) = ±
1. The complete figure set (3 figures)showing the same plot for cWB and LIB pipelines is availablein the online journal. ous parts of the sky. To quantify this feature, we intro-duce the extent ( δθ inj ) of a skymap as the maximum an-gular distance between the location of the injected signaland any other point satisfying p sky ( α, δ ) ≥ p . We showhistograms of δθ inj in Figure 3. The shown distributionsare clearly bimodal with peaks at cos( δθ inj ) = ±
1. Thepeak at cos( δθ inj ) = 1 corresponds to well localized sig-nals, while the peak at cos( δθ inj ) = − BBH sig-nals have twice as high peak at cos( δθ inj ) = 1 than thehistogram for G signals. Figure Set 3 shows histogramsof δθ inj obtained with the cWB and LIB pipelines on thesubset of signals identified as signals by BW.Even if previous measures indicate a well localizedsource, the skymap can still be fragmented, which makesit more difficult to cover the whole with EM observa-tions. We thus introduce the fragmentation of a skymapas the number of disjoint regions in the union of pointssatisfying p sky ( δ, α ) ≥ p . We show the distribution ofthe number of disjoint regions in Figure 4. Number ofdisjoint regions is less than 4 for more than 50% of in-jected signals for all morphology. Skymaps for SG and WNB signals are significantly more fragmented than for
Figure 4 . Distributions of fragmentation . Each row corre-sponds to one of the four morphologies ( SG , G , WNB , BBH ).Numbers at the bottom of the chart represent the number ofdisjoint regions in parts of the sky where p sky ≥ p . Numberof disjoint regions is less than 4 for more than 50% of in-jected signals for all morphologies. The complete figure set(3 figures) showing the same plot for cWB and LIB pipelinesis available in the online journal. G and BBH signals. This is due to the fact that theskymaps of these signals are more likely to have “fringepeaks”. These are separate rings in the sky correspond-ing to local maxima of matches between different datastreams obtained when they are shifted by half-integermultiples of the period of the signal (for details see Ap-pendix A). Figure Set 4 shows distributions of the num-ber of disjoint regions obtained with the cWB and LIBpipelines on the subset of signals identified as signals byBW.To compare BW’s performance with LIB’s and cWB’s(Essick et al. 2015), we created the equivalents of Fig-ures 1-4 with LIB and cWB using the same subset ofevents we used in this study (see Figure Sets 1-4). Wehave found that all metrics show that these algorithmsperform similarly in localizing the source. Histogramsof A show that A values for BW are comparable to,but systematically bigger than for cWB and LIB for allmorphologies, except for BBH signals, for which BWtypically yields smaller searched areas than LIB. Also,there are more
WNB skymaps with large searched ar-eas ( A &
100 deg ) for LIB than for BW. This is likelydue to its ability to recover more of the signal by usingmultiple wavelets as opposed to a single sine-Gaussiantemplate. 3.2. Waveform reconstruction
BW uses sine-Gaussian wavelets to reconstruct a GWsignal from the detector output, which means that therecovered signal is always given as a linear combinationof sine-Gaussian wavelets, the number of which is a pa-rameter in the RJMCMC. To characterize the qualityof waveform reconstruction, we introduce the overlap ( O , sometimes referred to as match) which measures Table 2 . Summary statistics of A and δθ distributions. Statistical errors are in the order of a few percent.morphology BBH SG G WNB fraction (in%) withsearched arealess than 5 deg . . . .
720 deg . . . . . . . . . . . . . . . . . . . . δθ less than 1 ◦ . . . . ◦ . . . . ◦ . . . . ◦ . . . . ◦ . . . . ◦ . . . . median δθ ◦ ◦ ◦ ◦ the similarity of an injected ( h i ) and a recovered ( h )waveform as: O = ( h i | h ) p ( h i | h i )( h | h ) , (5)where ( . | . ) is a noise weighted inner product, defined as:( a | b ) = 2 Z ∞ a ( f ) b ∗ ( f ) + a ∗ ( f ) b ( f ) S n ( f ) d f, (6)where S n is the one-sided power spectral density of thedetector noise, and x ∗ denotes the complex conjugate of x .From Eq. (5) it is visible that O ranges from -1 to 1,with O = 1 meaning perfect match between h i and h , O = 0 meaning no match at all, and O = − h i and h . With Eq.(5), we can calculate the overlap using data from onlyone detector. To characterize the waveform reconstruc-tion for the network of GW detectors, we introduce the network overlap ( O net ) by changing the inner productsin Eq. (5) with the sum of inner products calculated fordifferent detectors: O net = P Nj =1 ( h ( j )i | h ( j ) ) qP Nj =1 ( h ( j )i | h ( j )i ) · P Nj =1 ( h ( j ) | h ( j ) ) , (7)where j denotes the j -th detector in the network, and N is the number of detectors used in the analysis (notethat N = 2 in this study). Note that in our analysiswe only considered waveforms reconstructed from out-puts of each detector ( h ( j ) ), but not the astrophysicalGW polarizations ( h + , h × ), because the two polariza- tions cannot be decomposed from detections with twoco-aligned GW detectors, such as H1 and L1.Figure 5 shows the cumulative distribution functions(CDF) of O net . Shaded ranges represent the 2 σ uncer-tainty calculated using the Dvoretzky–Kiefer–Wolfowitzinequality (Dvoretzky et al. 1956). The fraction of in-jected signals with O net > . G , 96% for SG ,48% for BBH , and 47% for
WNB signals after the wave-form reconstruction with BW. 95% of injections have O net > .
92 for G signals, O net > .
91 for SG signals, O net > .
75 for
BBH signals, and O net > .
68 for
WNB signals. In Figure 5, the lower the curves reach at agiven O net value, the better the reconstruction is. Thissuggests that BW’s waveform reconstruction works mosteffectively for SG and G signals, for which the curves areidentical within the 2 σ statistical error. BW’s waveformreconstruction is less effective for WNB and
BBH sig-nals, and it shows similar characteristics for these mor-phologies at high network overlaps ( & . WNB signals has a longer tail at low O net values. The better performance of BW for SG and G signals is due to the fact that at low SNR net BW tends touse fewer wavelets to avoid overfitting the data. SG and G signals can be reconstructed accurately even with just2-3 sine-Gaussian wavelets, while this is not possible for WNB and
BBH signals. This also means that the curvesfor SG and G signals in Figure 5 represent BW’s maxi-mal capability of reconstructing a GW signal for a givennoise level, while the results for WNB and
BBH sig-nals represent BW’s performance on more generic (andthus, more realistic) GW signals. Note that while O net values are smaller for WNB and
BBH signals, BW de-tects these with more confidence, because its detection
Figure 5 . Cumulative distribution function (CDF) of net-work overlaps ( O net ). Shadings represent the 2 σ uncer-tainties calculated using the Dvoretzky–Kiefer–Wolfowitz in-equality (Dvoretzky et al. 1956). The lower the curves reachat a given O net value, the better the reconstruction is. Theinset shows the normalized histogram of network signal-to-noise ratio (SNR net ) for signals with four different morpholo-gies. The curves for SG and G signals are identical withinthe 2 σ statistical errors, and they indicate significantly bet-ter reconstructions of SG and G signals than of WNB and
BBH signals. statistic has a stronger dependence on signal complexitythan on SNR net (for details see Littenberg et al. 2016).The inset plot in Figure 5 shows the normalized his-togram of injected signals’ network signal-to-noise ratio(SNR net ) for the four different signal morphologies. SG and G signals have an overabundance at SNR net . WNB and
BBH signals. This indicates thatthe previously described difference in the distribution of O net is not due to the different SNR net distributions,as BW performs better for SG and G signals despitethe fact that SNR net values for SG and G signals areusually smaller than for WNB and
BBH signals. Notethat these distributions strongly depend on the parame-ter distributions of injected signals defined in Table 4of Essick et al. (2015), and on the corresponding de-tection efficiencies of the combination of cWB and BWpipelines for the different parameter sets (see the SNR net histogram in the inset of Figure 5).We show O net vs. SNR net for SG , G , and WNB sig-nals in the left panel of Figure 6. Curves were esti-mated with a Gaussian kernel smoother, which is a non-parametric regression method. Shaded regions betweendashed lines represent the 1 σ uncertainty regions calcu-lated with the bootstrap method, in which we estimatethe curve repeatedly for sub-samples randomly drawnfrom the full sample. Note that we excluded the injec-tions with SNR net >
100 from the estimation of thesecurves, and we only show the estimated curves up toSNR net =70. All three morphologies show a clear trend of O net increasing with SNR net .For BBH signals we calculated the O net vs. SNR net curves in two separate bins of total mass ( M tot ) of thebinary black hole system, calculated in the detector’sframe. The two bins were defined with M tot being M tot < ˆ M tot and M tot > ˆ M tot , where ˆ M tot = 44 . M ⊙ is the median of M tot values for all BBH injections. The O net vs. SNR net curves for BBH signals are shown inthe right panel of Figure 6. Similarly to other mor-phologies,
BBH injections also show a clear trend ofincreasing O net with increasing SNR net . At low ( . net values, BW performed significantly better forsignals with higher M tot , while differences in the curvesare within the level of statistical errors for higher SNR net values. Signals with high M tot are recovered with betteraccuracy because a large fraction of the signal power is ina compact region of time-frequency space and thereforecan be captured with a small number of wavelets, whilesignals with low M tot spend a comparatively longeramount of time in the sensitive band of the detectors, re-quiring more wavelets, and a larger total signal strength,to achieve a similar fit. This difference vanishes at highSNR net because BW uses more wavelets to reconstructsignals with higher SNR net .Similarly to Figure 5, Figure 6 also shows that BWperforms very similarly on SG and G signals, and muchless efficiently on WNB and
BBH signals. This is dueto the fact that BW needs to use more wavelets to accu-rately reconstruct
WNB and
BBH signals. Note, thatdespite the weaker performance on
WNB and
BBH sig-nals, these also approach the reconstruction accuracy for SG and G signals at higher SNR net values. Comparingthe two panels of Figure 6, it is visible that the curvefor BBH signals is similar to the curve for
WNB signals,with slightly worse overlap at low SNR net and slightlybetter overlap at high SNR net values.Our results show that BW robustly reconstructs wave-forms with various morphologies. Although there aresignificant differences between the efficiency of recon-structions of signals with different morphologies, evenfor the worst case of
WNB signals (which do not evenmatch BW’s assumption of the signal always being el-liptically polarized), most of them have relatively highoverlaps, and there is a clear trend of O net approaching1 as SNR net increases.3.3. Point estimates of waveform central moments
For a generic burst signal, we do not have any spe-cific astrophysical model whose parameters could be es-timated. In this case we can still give estimates onmodel-independent parameters of the signal. Here weconsider the central moments of the waveform as suchparameters.The first central moments are central time ( t ) and Figure 6 . Dependence of network overlaps ( O net ) on network signal-to-noise ratios (SNR net ) for SG , G , WNB , and
BBH signals. Note that we excluded the injections with SNR net >
100 from the curve estimation. Shaded areas represent the 1 σ uncertainty regions of the measured O net values. Left panel shows SNR net dependence of O net for SG , G , and WNB signals.All three morphologies show a clear trend of increasing overlap with increasing SNR net . Right panel shows SNR net dependenceof network overlaps for
BBH signals with detector-frame total mass below and above the median total mass ˆ M tot = 44 . M ⊙ .BW performed significantly better for signals with higher M tot at SNR net .
35 values. central frequency ( f ), and the second central momentsare duration (∆ t ) and bandwidth (∆ f ), defined as t = Z ∞−∞ d t ρ TD ( t ) t (8a) f = Z ∞ d f ρ FD ( f ) f (8b)(∆ t ) = Z ∞−∞ d t ρ TD ( t )( t − t ) (8c)(∆ f ) = Z ∞ d f ρ FD ( f )( f − f ) (8d)respectively, where ρ TD and ρ FD are the effective nor-malized distributions of signal energy, expressed in thetime domain (TD) and in the frequency domain (FD): ρ TD ( t ) = h ( t ) h , (9a) ρ FD ( f ) = 2( | ˜ h ( f ) | h , (9b)where h ( t ) is the whitened (i.e. normalized with the am-plitude spectral density of the detector noise) waveformfor a given detector and ˜ h ( f ) is the Fourier transform of h ( t ). These distributions satisfy R ∞−∞ ρ TD ( t )d t = 1, and R ∞ ρ FD ( f )d f = 1.Estimations of higher order moments could also begiven with BW, however we excluded them from our analysis due to the fact that they are more strongly af-fected by statistical errors than estimations of the firstorder moments (for a detailed discussion of this, see theend of this section).BW reconstructs the waveform and calculates thewaveform moments for each sample in the Markov chain.We calculated the median value to give a point estimateof the waveform moments. To quantify the accuracy ofthe point estimate of waveform moment x , we define theabsolute error of the estimation, e x , as: e x = | x (e) − x (r) | , (10)where x (e) is the estimated, and x (r) is the real value of x . We also introduce the relative error of an estimate, η x , as: η x = e x x (r) . (11)We show CDFs of e t / ∆ t , e f / ∆ f , η ∆ t , and η ∆ f inFigure 7, where shadings represent the 2 σ uncertain-ties calculated using the Dvoretzky–Kiefer–Wolfowitzinequality (Dvoretzky et al. 1956). All moments werecalculated for H1 detector data, however, results arevery similar for L1 too. We divided the absolute er-rors of the first moment estimations with the real valuesof the corresponding second moments, because we ex-pect that the statistical error of first moment estimationscales with the real values of the second moments.We show CDFs of e t / ∆ t for different morphologiesin the top left panel of Figure 7. These show that the0most accurate t estimates with BW are obtained for G signals, while estimates for BBH signals are the leastaccurate. The relatively large e t values are due to thefact that BW cannot reconstruct low-amplitude parts ofthe signal overwhelmed by noise, which can cause a sys-tematic error in the estimation of t . For example, BWis almost insensitive to the inspiral parts of BBH sig-nals, which make up the bulk of
BBH signal durations,and this bias increases the smaller the total mass of thesystems are. This effect is less significant for the otherthree morphologies, which explains why the estimationof t is less accurate for BBH signals (with a median e t value of 0 . t ). The t values we obtain for H1 and forL1 are strongly correlated, which means that the erroron the estimation of the difference of arrival times be-tween H1 and L1 (determining the thickness of the skylocalization triangulation ring) is typically smaller than e t .We show CDFs of η ∆ t in the top right panel of Figure7. These curves significantly differ for different wave-form morphologies. Regarding the median η ∆ t , ∆ t es-timation is the most accurate for SG signals (with amedian value of 0.06) and the least accurate for BBH signals (with a median value of 0.57). Note however,that median values do not tell about how long the tailsof η ∆ t distributions are. From the four morphologies,the η ∆ t distribution for the SG signals have the longesttail (see top right panel of Figure 7). For η ∆ t .
1, CDFvalues for
BBH signals are significantly smaller than forthe other three morphologies, while for η ∆ t &
1, theyare higher. This is due to the steep part of the
BBH curve around η ∆ t = 1, which corresponds to the system-atic underestimation of the duration of low mass BBH signals (due to the effect explained in the previous para-graph).We show CDFs of e f / ∆ f in the bottom left panel ofFigure 7. Curves for different morphologies are identicalwithin the error bars, in contrast with CDFs of e t / ∆ t where the curves are similar but not identical. This indi-cates that these errors are purely due to the statisticalerrors of central frequency estimation, determined bythe non-zero value of ∆ f . Note that all e f values aresmaller than ∆ f , and the median of e f is smaller than0.1 for all morphologies (see Table 3).We show CDFs of η ∆ f in the bottom right panel ofFigure 7. The accuracies of ∆ f estimation are simi-lar for different morphologies, but not as much as for η f / ∆ f . 95th percentiles are between 0.2 and 0.4 forthe different morphologies. Note that relative errors ofbandwidth estimations tend to be higher than of centralfrequency estimations. This is due to the fact that es-timations of second order moments inherit errors fromestimations of lower order moments (see Eq. (8c) and(8d)), and thus have higher statistical errors. We ex- pect that estimation of third and higher order momentswould have even bigger errors, and thus we restrict ourattention to examining only estimations of the first twomoments. Medians and 95th percentiles of errors foreach moment and for each morphology are shown in Ta-ble 3.As a summary, results presented in Figure 7 show thatthe distributions of errors for f and ∆ f are very similarfor different morphologies, while distributions of errorsfor t and ∆ t show significant differences between dif-ferent morphologies. This also means that while errorsof f and ∆ f estimations are purely statistical, errorsof t and ∆ t estimations include systematics as well.The latter is due to the fact that BW cannot recon-struct low-amplitude parts of a signal overwhelmed bynoise, which may result with a systematic error in theestimation of t and ∆ t . It is clear that the accuracyof moment estimation is affected by how accurately sig-nals are reconstructed. However, we see identical CDFsof e f / ∆ f for different morphologies, while these havedifferent O net distributions, which suggests that O net isnot a good indicator of BW’s moment estimation accu-racy. CONCLUSIONWe presented a comprehensive multi-aspect study onthe performance of BW, a Bayesian GW burst PEpipeline used by the LIGO-Virgo Collaboration for re-constructing GW burst signals and their parameters.We injected a large number of simulated signals withfour different morphologies (sine-Gaussians, Gaussians,white-noise bursts, and binary black hole signals) intosimulated O1 aLIGO noise to test BW’s performance inthree different aspects of PE: sky localization, waveformreconstruction, and estimation of waveform central mo-ments (for details on the methods we used see Section2).BW localizes sources with a level of accuracy com-parable for all four morphologies, with the median sep-aration of actual and estimated sky locations rangingfrom 25.1 ◦ to 30.3 ◦ (see Table 2), and median searchedarea ( A , see Eq. 4) ranging from 99.2 deg to 252.8deg (see Section 3.1). This is reasonable accuracy for atwo-detector network, and is comparable to accuraciesof other localization pipelines (cWB and LIB) studiedpreviously (Essick et al. 2015). Histograms of A (seeFigure 2) show that A values for BW are comparableto, but systematically bigger than for cWB and LIB forall morphologies. The exceptions are BBH signals, forwhich BW’s A values are systematically smaller. Notethat the runtime of cWB and LIB is much shorter thanof BW.BW reconstructs waveforms as a linear combinationof sine-Gaussian wavelets. To measure the goodness of1 Table 3 . Medians (50th percentiles) and 95thpercentiles of waveform central moment errorsfor the SG , G , WNB , and
BBH signal mor-phologies. P denotes the percentile rank of val-ues given in the corresponding table columns. P Signal morphology
SG G WNB BBH e t / ∆ t η ∆ t e f / ∆ f η ∆ f reconstruction, we used the network overlap ( O net , seeEq. 7), which quantifies the similarity between the in-jected and the reconstructed signals. We have foundthat BW reconstructs signals with O net > . G , 96% of SG , 45% of WNB , and 47% of
BBH sig-nals (see Section 3.2). We have also found that (seeFigure 6) O net increases rapidly with increasing SNR net ,reaching O net = 0 .
95 at SNR net ≈
14 for SG and G , atSNR net ≈
50 for
WNB , and at SNR net ≈
35 for
BBH signals. These results suggest that we can expect verygood reconstruction ( O net > .
95) for almost any sig-nal with high ( &
50) SNR net , and reasonably good re-construction ( O net > .
85) for almost any signal withmoderate ( &
20) SNR net .We also examined how accurately BW can estimatethe central moments of a GW waveform (see Section3.3). These are model-independent parameters of a sig-nal, therefore by examining the estimation of them, wecan characterize PE without assuming any astrophysicalmodel for the source. We have found that errors of f and ∆ f estimations are purely statistical, while errors of t and ∆ t estimations include some systematics as well.We have also found that O net is not a good indicator ofBW’s moment estimation accuracy. The median valueof e f / ∆ f is 0.09 for SG , G and WNB signals, and 0.07for
BBH signals (see Table 3). There is no standardprocedure on how the estimated moments of GW burstscan be used to test astrophysical models, however fu-ture studies can use our results to test the feasibility ofparticular methods using signal moments.This paper fits into a series of studies examin- ing PE for GW bursts (see e.g. Klimenko et al. 2011;Essick et al. 2015). These studies can be used incomparisons with improved performances of future PEpipelines, and in testing the feasibility of possible astro-physical applications of future GW burst detections.This paper was reviewed by the LIGO Scientific Col-laboration under LIGO Document P1600181. We thankMarco Drago and Sergey Klimenko for their valuablecomments on the manuscript. We acknowledge theBurst First2Years sky localization Open Data release .The authors acknowledge the support of the NationalScience Foundation and the LIGO Laboratory. LIGOwas constructed by the California Institute of Tech-nology and Massachusetts Institute of Technology withfunding from the National Science Foundation and op-erates under cooperative agreement PHY-0757058. Theauthors would like to acknowledge the use of the LIGOData Grid computer clusters for performing all the com-putation reported in the paper. Bence B´ecsy was sup-ported by the ´UNKP-16-2 New National Excellence Pro-gram of the Ministry of Human Capacities. Bence B´ecsywas supported by the Hungarian Templeton Programthat was made possible through the support of a grantfrom Templeton World Charity Foundation, Inc. Theopinions expressed in this publication are those of theauthors and do not necessarily reflect the views of Tem-pleton World Charity Foundation, Inc. Peter Raffai isgrateful for the support of the Hungarian Academy ofSciences through the ”Bolyai J´anos” Research Scholar-ship programme. Figure 7 . Cumulative distribution functions (CDF) of waveform central moment errors: absolute errors of central time estima-tions divided by signal durations ( e t / ∆ t , upper left), relative errors of duration estimations ( η ∆ t , upper right), absolute errorsof central frequency estimations divided by signal bandwidths ( η f / ∆ f , lower left), and relative errors of bandwidth estima-tions ( η ∆ f , lower right). Shadings represent the 2 σ uncertainties calculated using the Dvoretzky–Kiefer–Wolfowitz inequality(Dvoretzky et al. 1956). Colors indicate CDFs for signals with sine-Gaussian ( SG ), Gaussian ( G ), white noise burst ( WNB ),and binary black hole (
BBH ) morphologies. We give values of 95th percentiles and medians in Table 3.
APPENDIX A. RESOLVING BAYESWAVE O1 VERSION’S ISSUE WITH HIGH-Q SIGNALSDuring O1, BW was prone to classifying simulated short-duration high-frequency signals which underwent manywave cycles (i.e. high- Q signals) while in the measurement band of the detector as glitches. In principle there isno reason for the Bayesian evidence used to rank hypothesis under consideration by BW to have strong frequencydependence.Upon examination of the mis-classified injections, it was revealed that the high- f , high- Q signals exhibit multimodallikelihood support in the ( α, δ, t , f ) parameter sub-space. For these signals, the Markov Chain Monte Carlo (MCMC)sampler, which serves as the central engine to the BW algorithm, was not generically sampling between the differentmodes, and was thus prone to missing significant portions of the coherent signal and preferring the incoherent glitchmodel (which does not suffer the correlations between time-frequency parameters and sky location).3The cause of the multimodal likelihood function is clear. For a sinusoidal signal ( Q = ∞ ) the waveform is perfectlydegenerate when time-shifted by an integer number of wave-periods ( T ). For high- Q signals, a number of integerperiods (or half-integer periods with a π radians phase shift), time shifts produce similarly good fits to the data. Forcoherent signals, these (nearly) degenerate time shifts are also present in the time delay between detectors, which, forBW, is encoded in the sky location.To overcome BW’s susceptibility to missing modes of the likelihood when analyzing high- Q signals, we added aproposal distribution to the MCMC which explicitly suggests half-integer-period time shifts, along with half-integer-cycle phase shifts, for the wavelet parameters. Furthermore, extensive development (beyond the scope of this paper)to improve the overall capabilities of BW’s MCMC to sample the complicated sky-location posteriors encountered bytwo-detector gravitational-wave networks has been completed.
508 510 512 514 516 518 520 522 524 -6 -4 -2 0 2 4 6 f [ H z ] t-t [ms] φ [ r ad ] -1-0.5 0 0.5 1 0 1 2 3 4 5 6 s i n ( δ ) α [rad] -6-4-2 0 2 4 6 t - t [ m s ] Figure 8 . Scatter plot of MCMC samples for signal model parameters of a high- Q , high- f sine-Gaussian injection. The leftpanel shows the time-frequency plane with points colored by the wavelet phase parameter. Multiple modes and their phase-dependence are evident. The right panel shows the same chain samples, but now projected on the the sky-location plane ofthe parameter space and colored by the time parameter. Here again it is plain to see how different half-integer-period timeshifts correspond to different “rings” on the sky, making this a challenging distribution to sample without well-tuned proposaldistributions. Figure 8 contains two scatter plots from the BW MCMC utilizing the dedicated proposal distributions. The multi-modal nature of the posterior is on clear display, as is the efficiency with which the MCMC sampler is able to movebetween local maxima in the likelihood. This example came from an f ∼
512 Hz, Q ∼