Bayesian inference of overlapping gravitational wave signals
Elia Pizzati, Surabhi Sachdev, Anuradha Gupta, Bangalore Sathyaprakash
BBayesian inference of overlapping gravitational wave signals
Elia Pizzati,
1, 2
Surabhi Sachdev, Anuradha Gupta, and B. S. Sathyaprakash
1, 4, 5 Institute for Gravitation and the Cosmos, Department of Physics,Penn State University, University Park PA 16802, USA Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy Department of Physics and Astronomy, The University of Mississippi, Oxford MS 38677, USA Department of Astronomy and Astrophysics, Penn State University, University Park PA 16802, USA School of Physics and Astronomy, Cardi ff University, Cardi ff , CF24 3AA, United Kingdom (Dated: February 16, 2021)The observation of gravitational waves from LIGO and Virgo detectors inferred the mergers rates to be 23 . + . − . Gpc − yr − for binary black holes and 320 + − Gpc − yr − for binary neutron stars. These rates suggest that thereis a significant chance that two or more of these signals will overlap with each other during their lifetime in thesensitivity-band of future gravitational-wave detectors such as the Cosmic Explorer and Einstein Telescope. Thedetection pipelines provide the coalescence time of each signal with an accuracy O (10 ms). We show that usingthe information of the coalescence time, it is possible to correctly infer the properties of these overlapping signals with the current data-analysis infrastructure. Studying di ff erent configurations of the signals, we conclude thatthe inference is robust provided that the two signals are not coalescing within less than ∼ ∼ . ff er from significant biases in parameter inference, and newstrategies and algorithms are required to overcome such biases. I. INTRODUCTION
The advent of the third generation (3G) gravitational-wave(GW) observatories, such as the Cosmic Explorer (CE) [1, 2]and the Einstein Telescope (ET) [3], will o ff er the possibility toreach cosmological distances, thanks to an order of magnitudeimproved sensitivity compared to the current generation ofdetectors such as Advanced LIGO [4], Advanced Virgo [5],and KAGRA [6]. Indeed, 3G observatories will have unprece-dented sensitivity to detect coalescence events from an epochwhen the Universe was still in its infancy assembling its firststars and will routinely detect mergers with stupendously largesignal-to-noise ratios of several thousands [7, 8]. An order ofmagnitude greater redshift reach and access to extremely high-fidelity signals compared to current interferometers promisesmany new discoveries, while allowing completely independent,precision tests of cosmological models, alternative gravity the-ories, and astrophysical scenarios of compact binary formationand evolution [9]. With an expected rate of hundreds of thou-sands of binary coalescence signals each year on top of weakbut persistent radiation from isolated neutron stars, rare burstsfrom supernova and other transient sources and stochastic back-grounds, 3G observatories demand novel algorithms for signaldetection and characterization. Therefore, a proper understand-ing of systematics arising from overlapping loud and quietsignals alike will answer a range of scientific questions that areat the forefront of fundamental physics and astronomy, as wellas a realistic estimation of the computational cost.According to current estimates, 3G observatories are ex-pected to detect tens of thousands of binary black hole (BBH)and millions of binary neutron star (BNS) mergers each year[10, 11]. If we take account of the fact that signals will lastlonger due to a lower starting frequency (3 Hz for ET and 7 Hzfor CE), then it is clear that 3G data will be dominated bymany overlapping signals [12–15]. This poses two challenges:first, the detection of individual signals could, in principle, bea ff ected by the presence of multiple signals. Second, the cur-rent Bayesian inference methods [16, 17] may not guaranteeunbiased estimation of source parameters which is crucial to deliver the science promises of 3G observatories.The Laser Interferometer Space Antenna (LISA) is expectedto produce a data set containing many overlapping astrophysi-cal signals. Indeed, galactic white dwarf binaries are persistentsources of gravitational waves and they produce a “foreground”noise [18] that could masquerade the detection and parameterestimation of other astrophysical signals. Several authors havestudied the problem of both detection [19, 20] and Bayesianinference [21, 22] in this context, while others have focusedon the global solution to the full family of potential signals[23, 24] as well as implementation of novel algorithms that arespecifically suited to characterize multiple overlapping signals[25].The problem of overlapping signals producing a confusionbackground in future terrestrial detectors was identified morethan a decade ago [26]. Detecting overlapping GW signals hasbeen shown to be possible by two ET mock data challenges[12, 13]. These studies were able to correctly identify andrecover signals even when they were overlapping with multipleother signals. However, no e ff ort to study the problem ofinference in the case of terrestrial detectors has so far beenmade.In this work, we explore Bayesian parameter inference ofoverlapping signals in terrestrial detectors. Even though thesignal detection may provide unbiased results, there is noguarantee that the parameter inference in the case of over-lapping signals is possible within the current framework. Thisis because current methods heavily rely on the e ffi ciency ofsampling algorithms, which are used to explore the posteriordistribution of parameters. If we analyze overlapping signalswith the current parameter estimation (PE) procedures (i.e.,the assumption that the parameter space for multiple signalsis the same as in the case of data containing only one signalat a time), we expect Markov Chains and the posterior dis-tribution to exhibit a non-trivial behavior such as slowly ornon-convergence of chains, multi-modal and biased posteriordistributions, etc. To this end we will deploy the Fisher infor-mation matrix formalism, to gauge the limit between the regionwhere overlapping signals could lead to biases in parameter a r X i v : . [ g r- q c ] F e b inference and the region where they don’t. The Fisher studytells us that as long as the di ff erence in merger time ∆ t C oftwo overlapping signals is larger than the accuracy δ t withwhich their merger times can be measured (i.e., ∆ t C (cid:29) δ t ),irrespective of how long the individual signals are, parame-ter inference will not be a problem. We exploit this result inBayesian inference by choosing the prior on the merger epochas determined by the signal detection pipelines, which is about δ t C ∼ O (10 ms) . Indeed, most signals are recovered by searchpipelines with a bias of δ t C <
20 ms. A conservative prior onthe merger time could be a factor of 10 to 20 larger or at most0.5 ms. Thus, two overlapping signals with their merger timesseparated by larger than about 1–2 s do not su ff er from anysystematic biases. Hence, it su ffi ces to consider the extent towhich overlapping signals pose a problem for ∆ t C (cid:46) . The rest of the paper is organized as follows: in Sec. II,we compute the number of chunks in a year’s worth of datacontaining more than one merger. Section III is devoted tostudying the covariance between overlapping signals using theFisher information matrix with the emphasis on what we mightexpect for parameter inference in case of overlaps. Bayesianinference of overlapping signals is presented in Sec. IV. Ourmain conclusions and a brief discussion of the type of problemsthat should be addressed in future studies is presented in Sec. V.
II. NUMBER OF OVERLAPPING SIGNALS
The number of overlapping signals depends on (a) the typi-cal duration of signals and (b) the rate at which they arrive atthe detector. At the leading order, the length ξ of a coalesc-ing compact binary signal starting from a gravitational-wavefrequency f s until merger is given by ξ = (cid:16) G M / c (cid:17) − / ( π f s ) − / , (1)where G is Newton’s constant, c is the speed of light and thechirpmass M is related to the component masses m and m via M ≡ ( m m ) / / ( m + m ) / . A BNS system consistingof a pair of 1 . M (cid:12) would last for ξ (cid:39) s starting from afrequency of f s =
10 Hz (relevant for Advanced LIGO andAdvanced Virgo), 1.8 hr for f s = f s = M at acosmological redshift of z would appear in the detector to havea chirpmass of (1 + z ) M , and hence lives for a shorter durationin a detector’s sensitivity band. Thus, BNSs (1 M (cid:12) ≤ m , m ≤ M (cid:12) ) could last for tens of minutes to several hours in bandwhile BBH signals (3 M (cid:12) ≤ m , m ≤ M (cid:12) ) could last for tensof seconds to thousands of seconds.The cosmic merger rate of compact binary coalescencesdetermined by the first three observing runs of LIGO andVirgo [27, 28] implies that in a network of 3G observatoriesthe detection rate r , defined as the number of signals whosematched filter signal-to-noise ratio is larger than 12, lies inthe range r BBH ∈ [5 × , . × ] yr − for BBHs and r BNS ∈ [10 , ] yr − for BNSs [15]. Thus, given that signalslast for several hours, 3G data would contain several loudoverlapping signals. We shall see below that for the purposeof parameter inference the relevant quantity is not how manyoverlapping signals there are at any one time but if two or more signals have their merger times lie within a duration ∆ t . Thisis what we will set out to compute next.
A. Overlapping signals of the same family
Let r denote the Poisson detection rate of a given signalfamily (BBH or BNS). In an interval ∆ t , the expected Poissonrate is ν = r ∆ t and the probability of observing exactly k mergers during ∆ t is given by P k ( ν ) = ν k e − ν k ! . (2)Thus, the probability of observing two or more mergers during ∆ t is P k ≥ = ∞ (cid:88) k = P k ( ν ) = ∞ (cid:88) k = ν k e − ν k ! = − e − ν (1 + ν ) . (3)We have made use of the fact that the Poisson distribution isnormalized, namely (cid:80) ∞ k = P k ( ν ) = . To compute the numberof chunks N k ≥ in which two or more mergers will be observedwe must multiply the probability P k ≥ by the number of chunks n ∆ t = T / ∆ t in an observational period T : N k ≥ ≡ P k ≥ n ∆ t = (cid:2) − e − ν (1 + ν ) (cid:3) T ∆ t . (4)Substituting ∆ t = ν/ r and noting that N T ≡ r T is the totalnumber of signals detected during the period T , we get N k ≥ = (cid:2) − e − ν (1 + ν ) (cid:3) N T ν . (5)It is easy to see that in the limit ∆ t → ν → , N k ≥ (cid:39) ν N T / . The factor of 1 / ν. In the otherlimit, when ∆ t → T (and ν (cid:29) N k ≥ (cid:39) N k ≥ in which we canexpect to find two or more mergers in a year’s worth of data(i.e., using T = ν = r ∆ t ). Also indicated in the plot arethe detection rate of BBH (BNS) which is expected to be in therange r BBH ∈ [1 . , . × − s − ( r BNS ∈ [3 . , × − s − ,respectively) [15] in a 3G detector network comprising ofone ET and two CEs (one in north America and the other inAustralia). As we shall see in Sec. III, parameter inferenceshould not be a problem if the di ff erence in coalescence timesof a pair of signals is larger than ∼ ∆ t = . Thus, in Sec. IVwe will focus on Bayesian inference of signals whose mergertimes di ff er by about one second. We see that at the higherend of the BNS rate, we expect ∼ ,
000 one-second longchunks with two or more mergers while at the lower end ofthe BNS rate this number is ∼ . Likewise, ∼
300 chunkswill contain two or more BBH mergers at the higher end of theBBH detection rate while this number is ∼
40 at the lower endof the BBH rate.The detection rate of BBH signals in the current detectornetwork of LIGO, Virgo and KAGRA at their design sensitivity .
001 0 .
002 0 .
004 0 .
008 0 .
016 0 . r (in Hz)0 . . . . . . . . c hun k w i d t h t ( i n s ec o nd s ) t c hun k s w i t h m e r g e r s i n o n e y r
3G BNS detection rate .
3G BBH detection rate Bayesian inference of overlapping signals is unlikely to be a problem
Bayesian inference of overlapping signals could lead to systematic bias .
16 31600
FIG. 1. Contour diagram showing the number of times two or more signals have their epoch of coalescence occurring within an interval ∆ t in ayear’s worth of data as a function of the chunk size ∆ t and the Poisson rate r . Also shown are the detection rate of BBH and BNS signals in 3Gobservatories of one ET and 2 CEs [15]. As an example, if the detection rate is 8 mHz then we can expect in one year’s of data 1000 one-secondlong chunks in which two or more mergers would occur. For a pair of signals whose coalescence times di ff er by an interval of ∆ t > ∆ t < ∆ t → . is at best r ∼ . × − s − (or 730 yr − ) [28]. Thus, theprobability of observing multiple mergers in a chuck of size 1s or less is negligibly small in the Advanced detector era. Thiswill also be the case in the A + era [29] where the detectionrates are expected to be 3 times larger. B. Overlapping signals from two di ff erent families If the detection rate of signal families A and B are r A and r B , then probability that one or more mergers of each of thesesignal families would occur during an interval ∆ t is P A , k ≥ = − e − ∆ t r A , P B , k ≥ = − e − ∆ t r B . (6)Thus, the probability P AB that an interval ∆ t contains oneor more from each of the two signal families is simply theproduct P AB = P A , k ≥ P B , k ≥ . If the rates are small, this reducesto P AB = ( ∆ t ) r A r B and the number of such chunks over aperiod T is N AB = ( ∆ t ) r A r B T = N A N B / n ∆ t , where N A and N B are the total number of mergers during the period T offamilies A and B , respectively, and n ∆ t = T / ∆ t is the numberof chunks of width ∆ t during T . Using the range of BNS andBBH rates quoted before, we find that N AB would lie in therange 170–5100 for T = ∆ t = . From the foregoing discussions it is clear that a small butsignificant fraction of signals would have their coalescencetime within an interval of 1 s. As we shall see in the nextSection, due to their long duration, overlapping BNS signals are far less correlated with each other than overlapping BBHsignals. For the same reason, a pair of overlapping BNS andBBH signals are poorly correlated. Hence, in the Bayesianinference problem (Sec. IV) we will only consider overlappingBBH signals.
III. COVARIANCE AMONG OVERLAPPING SIGNALS
If two signals are well separated then the covariance betweentheir parameters is zero and we do not expect one signal toa ff ect the parameter inference of the other. As we bring thetwo signals closer together in time, at some point the presenceof one of the signals will begin to bias the estimation of param-eters of the other. In this Section we estimate the covariancebetween the parameters of a pair of overlapping signals usingthe Fisher matrix formalism. Although Fisher matrix is validin the limit of large signal-to-noise ratios, any inferences wecan draw from the correlation will guide us in choosing theparameter space of compact binaries where systematic biasescould be large.To this end, we assume that the data contains a pair of signals s A and s B buried in stationary, Gaussian noise n . The detectoroutput is a sum of the overlapping signals buried in detectornoise: x ( t ) = n ( t ) + s A ( t , λ ( A ) α ) + s B ( t , λ ( B ) α ) . (7)where λ ( A ) α , λ ( B ) α , for α = , . . . , p , are the set of parameters -2 -1 0 1-0.40.00.40.8 C o rr e l a ti on C o e ff i c i e n t σ M A M B σ η A η B σ t CA t CB -1 0 1 2 -0.4 0.0 0.4 0.8-2 -1 0 1Difference in merger times τ (seconds)-0.1-0.050.000.05 -1 0 1 2 -0.1-0.05 0.00 0.05BBH, Advanced LIGO BBH, Cosmic ExplorerBNS, Advanced LIGO BNS, Cosmic Explorer FIG. 2. Plot shows the correlation coe ffi cients, i.e., normalized covariances as defined by Eq. (16) between the parameters of the two overlappingsignals as a function of the di ff erence in merger times τ = t BC − t AC . The left panel is for Advanced LIGO and right for Cosmic Explorer. Toprow is for BBHs and bottom row BNS. Parameter inference of overlapping signals will be negligibly a ff ected when (the absolute value of) thecorrelation coe ffi cients are less than 10%. corresponding to signals s A and s B , respectively. Note thatsince both s A and s B are assumed to belong to the same sig-nal family they are specified by the same number of param-eters. Furthermore, we shall only consider a single detectorfor this exercise. The relevant parameters for a binary withnon-spinning companions are the chirpmass M , symmetricmass ratio η , the epoch t C when the signal amplitude reachesits peak and the phase φ C of the signal at that epoch and so: λ ( A ) α = ( M ( A ) , η ( A ) , t ( A ) C , φ ( A ) C ) and similarly for signal s B . Inthis Section, we assume the waveform model to be the T ay - lor F2 approximation at 3.5 post-Newtonian order [30–32].The epoch of merger is taken to be the time when the in-stantaneous gravitational-wave frequency reaches the value f ISCO = (6 / π M ) − , which corresponds to be the inner moststable circular orbit for Schwarzschild black holes.For the computation of the covariance matrix it is more con-venient to consider that the data contains only one signal, i.e.,the sum of the two signals s = s A + s B , and it is characterizedby a double number of parameters: θ a = λ ( A ) a for a = , . . . , p and θ a = λ ( B ) a − p for a = p + , . . . , p . For a noise backgroundthat is stationary and Gaussian the covariance matrix C , whichis inverse of the Fisher matrix Γ , is given by: C ab = Γ − ab , Γ ab = (cid:42) ∂ s ∂θ a , ∂ s ∂θ b (cid:43) . (8) Here the scalar product of two waveforms (or any pair offunctions of time for that matter) h and g is defined as (cid:104) h , g (cid:105) ≡ (cid:60) (cid:90) f high f low ˜ h ( f ) ˜ g ∗ ( f ) S h ( f ) d f , (9)where (cid:60) stands for the real part of the integral, ˜ h and ˜ g arethe Fourier transforms of the signals h and g, respectively, g ∗ denotes the complex conjugate of g and S h ( f ) is the one-sided noise spectral density of the detector. In our study wewill use either the noise spectral density of Advanced LIGO[4] or that of the Cosmic Explorer [9]. The lower frequencycuto ff f low is chosen to be 20 Hz for Advanced LIGO and5 Hz for Cosmic Explorer. The upper frequency cuto ff f high is assumed to be the larger of the inner-most stable circularorbit frequency of the two overlapping signals, i.e., f high = max[(6 / π M ) − , (6 / π M ) − ], where M and M are the totalmass of the two overlapping signals.The Fisher matrix contains interference terms of the follow-ing type: Γ α, β + p = (cid:42) ∂ s A ∂λ ( A ) α , ∂ s B ∂λ ( B ) β (cid:43) . (10)Covariances are of primary interest in this Section as they cantell us the degree to which the presence of one signal a ff ects theparameter inference of the other. In order to measure the extentof covariance we consider two sets of overlapping signals:1. overlapping BBHs with masses:( m ( A )1 , m ( A )2 ) = (21 M (cid:12) , M (cid:12) ) (11)( m ( B )1 , m ( B )2 ) = (33 M (cid:12) , M (cid:12) ) . (12)2. overlapping BNSs with companion masses:( m ( A )1 , m ( A )2 ) = (1 . M (cid:12) , . M (cid:12) ) (13)( m ( B )1 , m ( B )2 ) = (1 . M (cid:12) , . M (cid:12) ) . (14)Furthermore, in all cases we choose( t ( A ) C , φ ( A ) C ) = (0 , , ( t ( B ) C , φ ( B ) C ) = ( τ, π/ , (15)and vary τ over the range [ − ,
2] s . The covariances between the chirpmass, symmetric massratio and epoch of coalescence are plotted in Fig. 2 as a functionof the parameter τ for overlapping BBHs (top panels) andBNSs (bottom panels) for noise spectral densities of AdvancedLIGO (left panels) and Cosmic Explorer (right panels). Othercross-covariances are negligibly small and not shown. Whatwe plot are the normalized covariances, i.e., a combination ofthe correlation coe ffi cients defined as: σ ab ≡ C ab √ C aa C bb , a (cid:44) b . (16)This quantity is strictly bounded between − + . A correla-tion coe ffi cient of + − σ ab ∼ . ffi cients are negligibly smallbeyond | τ | ∼ . | τ | ∼ . Thecorrelation remains insignificant in the case of BNSs both inAdvanced LIGO and Cosmic Explorer except when τ (cid:39) . From this analysis, we conclude that the parameter inferenceof overlapping BNS signals is likely to be less severe than thatof overlapping BBH signals. We will therefore consider onlythe latter class of signals in the remainder of this paper.
IV. BAYESIAN INFERENCE OF OVERLAPPING SIGNALS
In this Section, we support the results we have derived usingthe Fisher information matrix formalism (Sec. III) with a fullBayesian inference procedure. With this parameter estima-tion (PE) process, we are able to fully explore the posteriordistribution of the parameters that generated the signals. Thisis important, because it allows us to confirm the presence (expected from the Fisher study) of distinct maxima in theposterior, one for each signal coalescing within the time chunkconsidered. Moreover, thanks to this numerical approach, wecan explore more carefully the region where biases are ex-pected (i.e., within ∆ t C < λ describing a compact binary coalescence (CBC) waveform h ( λ, t ), we can write the posterior distribution for λ as: P ( λ | x , h ) = π ( λ ) L ( x | λ, h ) Z ( x ) , (17)where x is the detector output. This posterior can be exploredby using a sampling algorithm (e.g., MCMC, nested sampling).As in Sec. III, assuming that the data x contains two overlap-ping signals s A ( signal A ) and s B ( signal B ), then it can bewritten as: x = n + s A + s B , (18)where n is the noise of the interferometer. Note that, in prin-ciple, to perform a Bayesian analysis of two or more over-lapping signals we should broaden the parameter space, e.g., θ = { λ A , λ B } , in order to account for the presence of multipleoverlapping signals. However, since running a sampling algo-rithm requires significant amount of computational resources,in most cases this is not required. In fact, as argued in Sec. III, ifthe signals’ coalescence times are wide apart we do not expectthe presence of one signals to influence posterior distributionof parameters of the other. For this reason, in what follows weconsider the parameter space of a single CBC signal. We willreturn on this point later on when discussing possible biasesarising because of this choice. A. Choice of signal families
As already mentioned, in this analysis we focus only onBBH signals. This choice is motivated by the fact that co-variances among overlapping BNS signals as demonstratedin Sec. III is negligibly small when their merger times di ff erby more than a few milliseconds—a duration during whichwe do not expect more than one merger to occur even in 3Gdetectors. Moreover, BNS signals last for several hours in 3Gdetectors and tens of minutes in Advanced LIGO and Virgoimplying Bayesian inference takes a formidable amount ofcomputational resources (although new algorithms are alreadyshowing the promise of greatly reducing the computationalrequirement [33]). Furthermore, we also restrict our analysisusing Advanced LIGO sensitivity. As argued before, LIGO isnot a ff ected by the problem of overlapping signals, because therate and the duration of the signals are far too small to createany overlap. Here we are not really interested in reproducing arealistic set of overlapping data, though; instead, we want tofocus on the parameter estimation process. To do so, there isno substantial advantage in using 3G mock data: we expectthat our conclusions will be valid even if they are based on theanalysis Advanced LIGO mock data.The parameters of the overlapping BBH signals used inBayesian inferences is the same as what we used in Sec. III: FIG. 3. Signals in the time domain, for four di ff erent values of the time shift τ . signal A , signal B , and the resulting overlapping signal areplotted. The waveforms are generated using the approximant imrphenomv
2. The two luminosity distances are fixed to d ( A ) L = , d ( B ) L = nonspinning BBHs with masses as given in Eq. (12) and coa-lescence times and phases as given in Eq. (15). We ignore theposition of the sources in the sky and their orientation relativeto the detectors (setting all angles to zero). We do, however,include in our analysis the luminosity distance d L of the source.The parameter space we use in our analysis is thus: λ = { m , m , φ C , t C , d L } Note that our choice of sky position is the worst case scenario,because we are considering the two sources to have the sameexact location in the celestial sphere. In reality, if overlappingsignals arrive from di ff erent directions in the sky, they willhave di ff erent phase coherence amongst a network of detectorsand thus easier to discriminate. Thus, since our choice of skyposition is the worst case scenario, the parameter estimationproblem can only be better when sky position and orientationare taken into account.To explore di ff erent configurations of the parameters, wevary the two luminosity distances of the sources d ( A ) L and d ( B ) L ,keeping the distance of one of the sources fixed to 1 Gpc andsetting the other at either 500 Mpc, 1 Gpc, or 2 Gpc. We alsovary the time shift τ , defined in Eq. (15) as the epoch coales-cence of signal B . The resulting variations in the parametersets are: τ = {− . , − . , − . , . , . } (19) d ( B ) L , d ( A ) L = {
500 Mpc , , } (20)With these choices, there are 25 di ff erent possible configu-rations, each of which is analyzed for Bayesian parameterinference.In the inference problem we use a signal model that accu-rately represents the BBH waveforms. To this end, we use the IMRP henomv ff to be 20 Hz,which is consistent with the minimum frequency used in theLIGO / Virgo PE. In Fig. 3, we plot the two waveforms in thetime domain, for the di ff erent configurations of the parame-ters. We vary only the time shift since changing the distancesof the signals results in a simple re-scaling of the amplitudesas we have neglected the e ff ect of redshift on masses. Theresulting overlapping waveform is plotted as well. In Table I,we compute the expected matched filter SNR for the di ff erentpossible configurations of the parameters, here we focus on thedistances, since the coalescence time does not a ff ect the SNRvalue. SNR d L = . d L = d L = ff erent values of the luminosity distances d L .Note that applying a time shift to the signals do not change the valueof the SNR. B. Setting up Bayesian inference runs
Having created the mock data with overlapping signals wenext focus on parameter inference. Our analysis uses twoLIGO interferometers, but our conclusions are not significantlya ff ected by this choice: considering a di ff erent detector networkwould simply result in di ff erent SNRs for the signals as we FIG. 4. Summary of the results for the set of 25 runs, each one with a di ff erent configuration of the parameters τ , d (A) L , and d (B) L . A runs areshown on the left panel, and
B runs are on the right panel. The recovered values for the masses m and m are plotted, together with theiruncertainties. The true values of the masses for signal A and signal B are highlighted with the dashed horizontal lines. Note that the points in theplots, referring to the same time shift τ , are slightly shifted with respect to their exact value of τ so that they do not overlap with each other. The τ = . are not focusing on the sky position of the source. Althoughthis could in principle change the heights of the peaks in theposterior distribution, we do expect it to influence their relativeratios significantly, and hence the PE process we consider isexpected to hold for any network.The data set consists of 4 s of mock data from the two LIGOinterferometers with a stationary, Gaussian noise background.4 s is large enough to span the full length of the longer signal.We use the package bilby [17, 34] to facilitate injection ofsignals into mock data that mimics the power spectral densityof Advanced LIGO.For Bayesian inference of the parameters of the injected sig-nals we use the bilby software package to perform parameterinference, running the dynesty sampler [35]. dynesty is a dy-namic, nested sampling algorithm [36, 37], which is well suitedfor our purposes because it quickly achieves convergence, butat the same time it is able to handle non-trivial, multi-modaldistributions better than MCMC-based algorithms. We allowthe sampler to explore the likelihood surface with respect toall the parameters except φ C , over which the likelihood is an-alytically marginalized, and d L , over which the likelihood isnumerically maximized. Marginalization over φ C and d L cor-rectly accounts for the e ff ects of the parameters φ C and d L onthe resulting 3-d posterior [38, 39]. Consequently, our finalposterior distribution is over the parameters m , m , and t C . C. Bayesian priors
At the beginning of the analysis, we have to set the priorson the various parameters. We consider a uniform prior onthe phase φ C , with periodic boundary conditions, a power-lawprior on the luminosity distance, p ( d L ) ∝ d α L with α = , anda uniform prior on the two masses m and m over the range[10 M (cid:12) ,
50 M (cid:12) ]. As for the coalescence time, selecting thebest possible prior turns out to be a game-changing strategy. Infact, running a simulation with a wide prior on the time t C thatspans the merger times of the two overlapping signals leads tosignificant problems: while one of the two signals is alwaysrecovered correctly, the other is completely ignored by thesampling algorithm. A wide prior on t C , therefore, would onlyallow us to infer the parameters of the louder signal, withoutaccess to the weaker one.However, as already pointed out, previous work suggests thatsignal can always be detected, even if they are overlapping, andtheir merger time correctly identified [12, 13]. The detectionof a signal allow us to know its epoch of coalescence with verylow uncertainty (at the order of 10 ms). For this reason, we canassume to know the time of coalescence of the two overlappingsignals with a good degree of accuracy, and constrain ourparameter space choosing a prior on the coalescence timewhich is centered on the (fiducial) true value of the time t C , with a width of 100 ms . In this way, for each of the signals wecan isolate the region of the parameter space where we expectto find the correct parameters. This choice allows us to recoverthe correct parameters for both signal A and signal B .Therefore, for each of the 25 injections, we run the Bayesianinference procedure two times: the first one (we refer to it as run A ) aims to recover the true values of the parameters of signal A ; to this end, since t (A) C = . Run B , on the otherhand, focuses on the signal B peak in the parameter space; thus,the prior is chosen to be centered in t C = τ . D. Results
In Fig. 4, we show the results for the recovered values ofthe parameters: results from run A are shown on the left panel,while run B ones are shown on the right panel. Note thatin the two panels only the values for the masses are shown:that is simply because, since we are analytically marginaliz-ing over φ C and d L , the only parameters for which we cancreate a posterior distribution are the two masses m and m ,and the coalescence time t C . However, the latter is alwaysrecovered correctly and with a very small error on its value(around 1 ms), as shown in the inset the right-hand panel ofFig. 4. The ‘violin’ plot was created from the distributionsof t C for a single time shift configuration, namely run B with τ = − . t C is determined tomilliseconds accuracy.Furthermore, we show the corner plots in Fig. 5 for tworandomly chosen runs ( run A on the right side and run B onthe left), created with the following choice of the parameters: τ = − . d (A) L = d (B) L = . The corner plots show thatthe results resemble the ones obtained when the data containsonly one signal: the posterior on the two masses presentsa degeneracy along the line of equal chirpmass, and the truevalues of the masses (and of the coalescence time) lie inside the2 σ credible region. Therefore, we confirm that – as expectedfrom discussion in Sec. III – using current PE algorithms it ispossible to treat the case of overlapping signals. We find thatthat setting the appropriate prior on the coalescence time t C , as determined by the search pipeline, is critical in determiningthe parameters of the signal without any bias.Looking at Fig. 4, we note that promising results are foundfor all the possible parameter configurations, except for thecase in which the time shift τ is equal to zero (i.e., the twosignals are coalescing at the exact same time). In fact, for τ (cid:44) . A runs (left panel), give results forthe parameters that are within one (sometimes two) sigma fromthe true values λ ( A ) = { m =
21 M (cid:12) , m =
15 M (cid:12) } . Similarly,in most of the considered cases, the B runs shown on the rightpanel give results for the parameters that are close to the truevalues λ ( B ) = { m =
33 M (cid:12) , m =
29 M (cid:12) } .If the time shift is zero (the corresponding runs are high-lighted with a grey shadowed box in Fig. 4), then both runsrecover only same signal. This means the the two-peaked pos-terior is not well sampled by the algorithm. The chains are ableto find only one peak, and they are not able to sample e ffi cientlythe other one: a possible explanation for this behaviour is thatthe two peaks have very di ff erent values (on average, the twopeaks di ff er for 100 −
500 orders of magnitude), thus samplingthe minor peak turns out to be impossible for the algorithm. This hypothesis is supported by an examination of Table I: thehighest peak is the one that’s associated with the highest SNR.As expected, then, in the runs in which signal B is louder than signal A , we recover a peak close to λ ( B ) , and vice-versa.We note here that the fact that the algorithm was not ableto focus on both peaks at the same time was already clearfrom the results in the case of a wide prior on the coalescencetime: if we do not isolate the correct value of t C by using anarrow t C prior, again we recover only the loudest peak andits corresponding time shift, with no sign of the other peak inthe posterior. We could also try to isolate one peak at a timeby imposing a narrower prior on the two masses m and m .This is more di ffi cult though, because of the relatively largeuncertainties of the recovered values of the masses (see e.g.,Fig. 5). For this reason, it is not easy to isolate one peak bypreventing any influence of the other one. Moreover, imposinga narrower prior on the component masses does not have aclear justification as signal detection pipelines return fiducialvalues for the two masses only within a large credibility region,and this could prevent in principle a reliable narrowing of themass priors.In some particular cases, the recovered masses are morethan two-sigma away from the expected values of m and m . This behavior is more common for τ = . τ = . τ . In order to assess the reasons for this wrong behavior,we run the biased configurations without the presence of theinterferometer noise. This allows us to discern between theinfluence of the noise from the one linked to the presenceof an overlap in the data. We find that, in the case of theruns with τ (cid:44) . τ = . τ = . V. DISCUSSION AND OUTLOOK
We presented a Bayesian inference analysis in the case ofoverlapping gravitational waves signals. Our goal was to assessthe capabilities of current Bayesian inference infrastructure tohandle the non-trivial case of one or multiple overlaps happen-ing within a data segment. This problem is destined to play amajor role in 3G detector planning, since the dramatic increasein sensitivity will result in a great number of signals coalescingwithin a few seconds.We started from a study based on the Fisher matrix for-malism, in which we analyzed the correlation between two
FIG. 5. Corner plots for a single run A (left side) and a single run B (right side); the overlapping signals are created with the following choice ofthe parameters: τ = − . d (A) L = d (B) L = φ C and the luminosity distance d L , sothat the three parameters considered here are the two masses m and m , and the coalescence time t C . The true values of these parameters arehighlighted with red dashed lined in the corner plots. The blue dashed vertical lines represent the 1 σ error on the parameters. overlapping signals. In this way, we were able to determinewhether in some regions of parameter space the two signalsare strongly correlated between each other, thus preventinga distinct inference procedure for one signal at a time. Wefound that BNS signals are not expected to be significantlycorrelated, and therefore their inference could be a problemonly for coalescence times really close between each other (atthe 10 −
100 ms level). BBHs, instead, su ff er from the presenceof a correlation starting from a much greater time shift τ (i.e.,the di ff erence between the two coalescence times). However,also for BBHs we find that for time shifts wider than ∼ dynesty sampling algorithm to describe the posterior distribution forthe parameters considered. We showed that, in order to samplea single peak without worrying for the presence of the otherone, a possible solution is to impose a narrow prior around thefiducial value (provided by the signal detection pipeline) ofthe coalescence time of the signal of interest. This procedureallows to isolate one single peak at a time, and worked wellin the configurations we explored. However, as the time shiftapproaches zero, we witnessed the emergence of biases thatprevented the recovery of the true values of the parameters.These biases happened only in the zero time-shift case (al-though our time resolution was at the 0 . . < | τ | < . ff erent approach that wedid not attempt in this work. One possible solution is to simplybroaden the parameter space searching for multiple signals inthe same Bayesian inference run. This significantly increasesthe computational costs of the Bayesian algorithms; however, since such an algorithm would be needed only in the worst casescenario of | τ | ∼
0, then the number of overlaps in such a smalltime region will be rarely greater than two. This means that,in the light of our work, heavy Bayesian inference methodswith a very large number of signals (and parameters) are notreally necessary in the case of overlapping signals. Anotherpossible solution to the biases would be to create an iterativeprocedure where one hierarchically determines the parametersof louder signals (as inferred from search algorithms) andsubtracts them from the data before analysing weaker ones[11, 40]. Although it is not possible to use narrower priorson component masses to isolate multiple peaks, it should bepossible to impose constraints on the chirpmass as found bydetection pipelines. It is, however, important to ascertain theextent to which such constraints can imposed by carrying outthe detection problem on a large sample of injections and theaccuracy with which detection pipelines are able to measurechirpmass. Our study did not include other source parameterssuch as companion spins, the position of the source in the skyand the orientation of binary relative to the detector frame.In particular, when overlapping signals arrive from di ff erentpositions in the sky then they would, in general, have di ff erentcoalescence times in di ff erent detectors, which might help toisolate one of the peaks better [41] but we need to study thisproblem more carefully. These and related problems will beexplored in a future study. ACKNOWLEDGMENTS
We thank Anuradha Samajdar, Justin Janquart, Chris VanDen Broeck and Tim Dietrich for sharing and discussing theirresults on a similar study [15]. We are indebted to useful com-ments by Rory Smith and Salvatore Vitale. EP is grateful to0Walter Del Pozzo for useful suggestions. EP was supportedby INFN in the framework of the 2019 NSF-INFN Summerexchange program. SS is supported by the Eberly PostdoctoralFellowship of Penn State. BSS is supported in part by NSF Grant No. PHY-1836779. The authors are grateful for compu-tational resources provided by the LIGO-Caltech ComputingCluster. This paper has the LIGO document number P2100044. [1] LIGO S cientific collaboration,
Exploring the Sensitivity of NextGeneration Gravitational Wave Detectors , Class. Quant. Grav. (2017) 044001 [ ].[2] D. Reitze et al., Cosmic Explorer: The U.S. Contribution toGravitational-Wave Astronomy beyond LIGO , [ ].[3] M. Punturo et al.,
The Einstein Telescope: A third-generationgravitational wave observatory , Class. Quant. Grav. (2010)194002.[4] LIGO S cientific collaboration, Advanced LIGO , Class. Quant.Grav. (2015) 074001 [ ].[5] VIRGO collaboration, Advanced Virgo: a second-generationinterferometric gravitational wave detector , Class. Quant. Grav. (2015) 024001 [ ].[6] KAGRA collaboration, KAGRA: 2.5 Generation InterferometricGravitational Wave Detector , Nat. Astron. (2019) 35[ ].[7] B. Sathyaprakash et al., Scientific Objectives of EinsteinTelescope , Class. Quant. Grav. (2012) 124013 [ ].[8] S. Vitale and M. Evans, Parameter estimation for binary blackholes with networks of third generation gravitational-wavedetectors , Phys. Rev.
D95 (2017) 064052 [ ].[9] D. Reitze et al.,
The US Program in Ground-BasedGravitational Wave Science: Contribution from the LIGOLaboratory , Bull. Am. Astron. Soc. (2019) 141[ ].[10] V. Baibhav, E. Berti, D. Gerosa, M. Mapelli, N. Giacobbo,Y. Bou ff anais et al., Gravitational-wave detection rates forcompact binaries formed in isolation: LIGO / Virgo O3 andbeyond , Phys. Rev.
D100 (2019) 064060 [ ].[11] S. Sachdev, T. Regimbau and B. Sathyaprakash,
Subtractingcompact binary foreground sources to reveal primordialgravitational-wave backgrounds , Physical Review D (2020).[12] T. Regimbau, T. Dent, W. Del Pozzo, S. Giampanis, T. G. Li,C. Robinson et al.,
Mock data challenge for the einsteingravitational-wave telescope , Physical Review D (2012)122001.[13] D. Meacher, K. Cannon, C. Hanna, T. Regimbau andB. Sathyaprakash, Second einstein telescope mock data andscience challenge: Low frequency binary neutron star dataanalysis , Physical Review D (2016) .[14] T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis,B. Sathyaprakash and S. Vitale, Digging deeper: Observingprimordial gravitational waves below the binary black holeproduced stochastic background , Phys. Rev. Lett. (2017)151105 [ ].[15] A. Samajdar, J. Janquart, C. V. D. Broeck and T. Dietrich,
Biases in parameter estimation from overlapping gravitationalwave signals in the third generation detector era , [ ].[16] J. Veitch et al.,
Parameter estimation for compact binaries withground-based gravitational-wave observations using theLALInference software library , Phys. Rev.
D91 (2015) 042003[ ].[17] G. Ashton, M. H¨ubner, P. D. Lasky, C. Talbot, K. Ackley,S. Biscoveanu et al.,
Bilby: A user-friendly bayesian inferencelibrary for gravitational-wave astronomy , The AstrophysicalJournal Supplement Series (2019) 27. [18] J. Crowder and N. J. Cornish,
LISA source confusion , Phys. Rev.D (2004) 082004 [ gr-qc/0404129 ].[19] N. J. Cornish and E. K. Porter, The Search for supermassiveblack hole binaries with LISA , Class. Quant. Grav. (2007)5729 [ gr-qc/0612091 ].[20] T. B. Littenberg, A detection pipeline for galactic binaries inLISA data , Phys. Rev. D (2011) 063009 [ ].[21] N. J. Cornish and J. Crowder, LISA data analysis using MCMCmethods , Phys. Rev. D (2005) 043005 [ gr-qc/0506059 ].[22] J. Crowder and N. Cornish, A Solution to the GalacticForeground Problem for LISA , Phys. Rev. D (2007) 043008[ astro-ph/0611546 ].[23] T. Littenberg, N. Cornish, K. Lackeos and T. Robson, GlobalAnalysis of the Gravitational Wave Signal from GalacticBinaries , Phys. Rev. D (2020) 123021 [ ].[24] T. Robson and N. Cornish,
Impact of galactic foregroundcharacterization on a global analysis for the LISA gravitationalwave observatory , Class. Quant. Grav. (2017) 244002[ ].[25] A. Petiteau, S. Babak, A. Sesana and M. de Ara´ujo, Resolvingmultiple supermassive black hole binaries with pulsar timingarrays II: genetic algorithm implementation , Phys. Rev. D (2013) 064036 [ ].[26] T. Regimbau and S. A. Hughes, Gravitational-wave confusionbackground from cosmological compact binaries: Implicationsfor future terrestrial detectors , Phys. Rev. D (2009) 062002[ ].[27] LIGO S cientific , V irgo collaboration, Binary Black HolePopulation Properties Inferred from the First and SecondObserving Runs of Advanced LIGO and Advanced Virgo , Astrophys. J. (2019) L24 [ ].[28] LIGO S cientific , V irgo collaboration,
Population Properties ofCompact Objects from the Second LIGO-VirgoGravitational-Wave Transient Catalog , [ ].[29] KAGRA, LIGO S cientific , VIRGO collaboration,
Prospects forObserving and Localizing Gravitational-Wave Transients withAdvanced LIGO, Advanced Virgo and KAGRA , Living Rev. Rel. (2018) 3 [ ].[30] B. S. Sathyaprakash and S. V. Dhurandhar, Choice of filters forthe detection of gravitational waves from coalescing binaries , Phys. Rev. D (1991) 3819.[31] T. Damour, B. R. Iyer and B. S. Sathyaprakash, A Comparisonof search templates for gravitational waves from binary inspiral , Phys. Rev. D (2001) 044023 [ gr-qc/0010009 ].[32] A. Buonanno, B. Iyer, E. Ochsner, Y. Pan and B. Sathyaprakash, Comparison of post-Newtonian templates for compact binaryinspiral signals in gravitational-wave detectors , Phys. Rev. D (2009) 084043 [ ].[33] D. Finstad and D. A. Brown, Fast Parameter Estimation ofBinary Mergers for Multimessenger Follow-up , Astrophys. J.Lett. (2020) L9 [ ].[34] I. Romero-Shaw et al.,
Bayesian inference for compact binarycoalescences with BILBY: Validation and application to the firstLIGO–Virgo gravitational-wave transient catalogue ,[ ].[35] J. S. Speagle, dynesty: a dynamic nested sampling package forestimating bayesian posteriors and evidences , Monthly Notices of the Royal Astronomical Society (2020) 3132–3158.[36] J. Skilling, Nested Sampling , .[37] E. Higson, W. Handley, M. Hobson and A. Lasenby,
Dynamicnested sampling: an improved algorithm for parameterestimation and evidence calculation , Statistics and Computing (2018) 891–913.[38] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Gra ff , S. Vitale et al., Parameter estimation for compact binaries with ground-basedgravitational-wave observations using the lalinference software library , Phys. Rev. D (2015) 042003.[39] L. P. Singer and L. R. Price, Rapid Bayesian positionreconstruction for gravitational-wave transients , Phys. Rev. D (2016) 024013 [ ].[40] C. Cutler and J. Harms, BBO and the neutron-star-binarysubtraction problem , Phys. Rev. D (2006) 042001[ gr-qc/0511092 ].[41] P. Christian, S. Vitale and A. Loeb, Detecting Stellar Lensing ofGravitational Waves with Ground-Based Observatories , Phys.Rev. D (2018) 103022 [1802.02586