How to search for multiple messengers -- a general framework beyond two messengers
HHow to search for multiple messengers - a general framework beyond two messengers
Do˘ga Veske, ∗ Zsuzsa M´arka, Imre Bartos, and Szabolcs M´arka Department of Physics, Columbia University in the City of New York, New York, NY 10027, USA Columbia Astrophysics Laboratory, Columbia University in the City of New York, New York, NY 10027, USA Department of Physics, University of Florida, PO Box 118440, Gainesville, FL 32611-8440, USA
Quantification of the significance of a candidate multi-messenger detection of cosmic events isan emerging need in the astrophysics and astronomy communities. In this paper we show that amodel-independent optimal search does not exist, and we present a general Bayesian method for theoptimal model-dependent search, which is scalable to any number and any kind of messengers, andapplicable to any model. In the end, we demonstrate it through an example for a joint gravitationalwave, high-energy neutrino, gamma-ray burst event search; which has not been examined heretofore.
I. INTRODUCTION
Astronomy has started via observations made in thevisible region of the electromagnetic spectrum in ancienttimes [1, 2]. As the technology and physics knowledgeof humanity developed, more and better observationswere made with new equipment and via new messengers;such as the whole electromagnetic spectrum [3–8], cos-mic rays[9, 10], neutrinos [11] and recently gravitationalwaves [12]. The new messengers have made it possibleto observe events which had not been possible before,as well as to gather a more complete picture of a singleevent by probing different processes of it. This allows usto understand the ongoing physics at extreme conditionsthat we cannot produce on Earth.Three observations, each involving at least two mes-sengers, can be given as examples for multi-messengerdiscoveries. The first one was the supernova SN1987 ob-served in electromagnetic waves and low-energy neutri-nos (in MeV energy range) in 1987 [13]. The secondwas the observation of the binary neutron star merger,GW170817, which was first discovered with gravitationalwaves and gamma-rays [14]. Later it was tracked in all ofthe electromagnetic spectrum. Finally, the last one wasa flaring blazar observed in gamma rays and high-energyneutrinos[15].As detectors improve for all messengers, it is naturalto expect to have more multi-messenger detections withmore messengers and better data. Therefore a need fora framework for multi-messenger coincidence quantifica-tion is inevitable. For example, the HAWC observatoryrecently observed a subthreshold gamma-ray candidatecoming from the coincident sky area of a neutrino de-tected by IceCube in response to a significant gravita-tional wave detection by LIGO and Virgo detectors [16].One challenge here is relating different messengers ofthe same source to each other. The possibility of havingseveral unrelated detections or noise triggers coinciden-tally showing up in the appropriate spatial and tempo-ral regions for a potential multi-messenger observation ∗ [email protected] makes it impossible to deduce the multi-messenger de-tection with absolute certainty. Therefore a statisticalinference has to be made [17–21].In this paper, we first describe the main challenge for amulti-messenger search in Sec. II and show that a model-independent optimal search does not exist. We provide aBayesian solution for assigning a significance to a multi-messenger detection, or candidate observations of differ-ent messengers in Sec. III, which is the extension of themethod described in [17] for coincident high-energy neu-trinos and gravitational waves. In Sec. IV we demon-strate the method for a joint gravitational wave, high-energy neutrino, gamma-ray burst event search, whichhas not been examined until this work. We note that thedescribed method is scalable to any number or any typeof messengers. We conclude in Sec. V. II. THE MULTI-MESSENGER SEARCHPROBLEM
The problem we want to address in this paper is toconstruct an optimal search for multi-messenger events.These searches can be described as the analyses thatquantify the chance of a number of messengers comingfrom the same source. As that number can be at leasttwo (i.e., what is the chance that at least two of themessengers have come from the same source?), one canlook for multi-messenger events with a different numberof messengers, and can put a constraint to the type ofthe messengers as well.In terms of statistics, the problem for these searchesis a composite hypotheses testing problem. Our inputparameters are the detection properties of the messen-gers, which may or may not be of astrophysical origin.Correspondingly, let’s consider a search with n ( n ≥ m (0 ≤ m ≤ n ) of the messengers beingastrophysical and coming from l (1 ≤ l ≤ m ) existingastrophysical sources. Naturally, all of them being noiseoriginated ( m = l = 0) is also a possibility. The totalnumber of sub-hypotheses for n messengers is given by a r X i v : . [ a s t r o - ph . H E ] O c t f ( n +1), for the f function defined recursively in Eq. (1). f ( n ) = n − (cid:88) i =0 (cid:18) n − i (cid:19) f ( i ) , f (0) = 1 (1)For example, for two messengers, there are f (3) = 5sub-hypotheses, which are; both of them being not real(noise), only the first one being real, only the secondone being real, both of them being real and coming fromthe same source, and finally both of them being real andcoming from different sources.In the context of the multi-messenger search, the pos-sible sub-hypotheses form two distinct hypotheses, com-monly named null and alternative hypotheses. We willcall our alternative hypothesis as the signal hypothesis.For a multi-messenger search for at least two messengerscoming from the same source, the null hypothesis con-sists of the sub-hypotheses which have l = m , so thatnone of the astrophysical messengers have come from thesame source. The signal hypothesis contains the remain-ing sub-hypotheses. As there is a composite hypothesestesting problem, one may naturally look for the uniformlymost powerful (UMP) test, which does not exist for ourproblem in general as it will be illustrated. The UMP testis the level α test (false alarm or type I error probabilityfor all null sub-hypotheses is at most α ) which has thehighest statistical power (the least false dismissal or typeII error probability) for all of the signal sub-hypotheses.If we show that the level α most powerful tests for anytwo signal sub-hypotheses are different, then we will haveproved that the UMP test does not exist. It should benoted that search for UMP tests is meaningful when wehave more than two signal sub-hypotheses, as for a sin-gle signal sub-hypothesis one can always find the mostpowerful test. Hence we construct our illustration for n >
2, for having more than one signal sub-hypothesis.Consider the search for at least two messengers from thesame source with n = 3 with messengers: M , M , and M . The most powerful test for the signal sub-hypothesiswhich has M and M coming from the same source and M being unrelated to them favors the events which havespatial overlap between the localization of M and M ,for example the test statistic which is proportional to theproduct of M and M ’s 2D sky (or 3D volume if avail-able) localizations without involving M ’s localization inthe integral. However the most powerful test for the cor-responding signal sub-hypothesis for M and M comingfrom the same source and M being unrelated favors theevents which have spatial overlap between the localiza-tion of M and M . So the two most powerful tests can-not be same and UMP test for this search doesn’t exist.One advantage of dealing with the joint observationsof previously individually studied messengers is know-ing both their astrophysical and noise originated rate ofoccurrences. By using these rates, one can empiricallyweight and combine the sub-hypotheses which has l = m and also the sub-hypotheses involving multiple same typeof messengers coming from the same source. This re- duces the number of sub-hypotheses from f ( n + 1) to f ( n ). Combining more sub-hypotheses requires assum-ing rates for multi-messenger observations. Due to thesmall number of such detections, these rates could notbe empirically determined and cannot be used in an ob-jective manner. III. BAYESIAN STRATEGIES FOR MODELS
As discussed in the previous section, one has to make amodel-dependent choice to combine the remaining f ( n )sub-hypotheses after using the individual messenger de-tection rates, which provides the most powerful searchfor the chosen model. The ratio of the predicted num-ber density of multi-messenger sources to the numberdensity of individual messengers’ sources together withsources’ emission models (i.e. emission energies, depen-dency on inclination etc.) and messengers’ propagationsin space, the ratios between the rates of each kind of de-tection can be found, which is necessary for weighting allof f ( n + 1) sub-hypotheses. After weighting, null andsignal hypotheses reduce to simple hypotheses and theNeyman-Pearson lemma [22] can be used for finding themost powerful test. The resulting test statistic (TS) isgiven in Eq. (2).TS( x ) = P ( x | H s ) P ( x | H n ) = (cid:80) i P ( x | H is ) P ( H is ) (cid:80) j P ( x | H jn ) P ( H jn ) × (cid:80) j P ( H jn ) (cid:80) i P ( H is )(2)where x is the complete set of detection outcomes, H s and H n are the signal and null hypotheses in order,and H is and H jn are the individual signal and null sub-hypotheses in order. We will ignore the very last term inthe Eq. (2) as it does not depend on x . The hypothe-sis prior probabilities P ( H ) will be canceled with a sameterm in detection likelihoods P ( x | H ba ). A. Detection likelihoods P (x | H ba ) Next, we explain the detection likelihoods. There aretwo parts to the issue. The first one is the decision ofthe origin, whether a messenger is astrophysical or noiseoriginated. This part is generally decoupled for differ-ent types of messengers due to independent detectors.The detection outcomes for each messenger are used to-gether with the detector characteristics to determine thispart. The second one is the multi-messenger aspect of thedetection for which correlations between messengers arerequired, especially in the space-time coordinates of themessengers. This coupling can be done with a sourcemodel with parameters θ as P ( x | H ba ) = (cid:90) P ( x | θ , H ba ) P ( θ | H ba ) d θ (3)The source parameters θ can include any property ofthe sources (there can be more than one source depend-ing on the sub-hypothesis) such as emission energies orspatial position of the sources. Prior information of suchproperties can be summarised in a joint density distribu-tion P ( θ ). If the corresponding sub-hypothesis H ba doesnot include a multi-messenger detection, then there maynot be a requirement for a common source and the sourceparameters θ . In that case, if the detectors are indepen-dent from each other, the detection outcomes’ probabili-ties can be expanded as a product. P ( x | H ba ) = (cid:89) i P ( x i | H ba ) (4)where subscript i runs over different detectors and x i are the detection outcomes from i th detector. Similarly,when there is a common source we can expand the de-tection outcomes’ probabilities for a fixed source as aproduct for different detectors. P ( x |{ θ } , H ba ) = (cid:89) i P ( x i |{ θ } , H ba ) (5)There we used the notation { θ } = { θ , θ , ... } for rep-resenting possible set of sources. There can be an addi-tional level of complication related to combinatorics if thesubhypothesis H ba can be satisfied with different group-ings of the detections. To illustrate this, consider thesubhypothesis of having five detected particles, three ofthem from a source and all the rest being noise origi-nated. In this case the subhypothesis can be satisfiedwith (cid:0) (cid:1) = 10 combinations. In such cases we expandthe probabilities P ( x i | θ , H ba ) as P ( x i |{ θ } , H ba )= (cid:88) { x ji , x ki ,... } P ( x i |{ θ } , H ba , {{ x ji , x ki , ... } , { x pi , x qi , ... } , ... } ) × P ( {{ x ji , x ki , ... } , { x pi , x qi , ... } , ... }|{ θ } , H ba ) (6)where the sum is over all the combinations of de-tection outcomes satisfying the sub-hypothesis H ba ,the sets { x ji , x ki , ... } and { x pi , x qi , ... } are the detec-tion outcomes of individual detections from differentsources and P ( {{ x ji , x ki , ... } , { x pi , x qi , ... } , ... }|{ θ } , H ba ) isequal to the reciprocal of the total possible com-binations for H ba arising from detector i . For ex-ample, for a subhypothesis with a single sourceand w particles emitted from that source, and to-tal of W detections; P ( { x i , x i ... x wi }| θ , H ws ) = (cid:0) Ww (cid:1) − . P ( x i |{ θ } , H ba , {{ x ji , x ki , ... } , { x pi , x qi , ... } , ... } ) are found byphysics and empirical data. For memoryless detectors,the detections from different sources are independent so P ( x i |{ θ } , H ba , {{ x ji , x ki , ... } , { x pi , x qi , ... } , ... } )= P ( { x ji , x ki , ... }| θ , H ba ) P ( { x pi , x qi , ... }| θ , H ba ) (7)The likelihoods P ( { x ji , x ki , ... }| θ , H ba ) are found via thedetector characteristics and emission models. Now we look at our second term in Eq. (3), P ( θ | H ba ).We transform it by using the Bayes’ rule. P ( θ | H ba ) = P ( H ba | θ ) P ( θ ) P ( H ba ) (8)The denominator of Eq. (8) cancel with the same termin Eq. (2). P ( θ ) is the joint density of source parametersbeing integrated over.The sub-hypothesis probabilities P ( H ba | θ ) are foundvia the expected counts from the sources or the noiseorigin by assuming a source density and using the empir-ically known noise trigger rate. IV. USE CASES – EXAMPLE: A JOINTGRAVITATIONAL WAVE – HIGH-ENERGYNEUTRINO – GAMMA-RAY BURST EVENTSEARCH
The method for multi-messenger searches introducedabove can be used in all scenarios. Specifically, in high-energy astrophysics, one can search for sources whichemit more than one messengers. Those messengers canbe in any wavelength in the electromagnetic spectrum,can be neutrinos, cosmic rays, gravitational waves, or anyother messenger. Joint emissions of gravitational waves,high-energy neutrinos, and gamma rays from a binaryneutron star merger or a binary black hole merger in adense medium, such as an AGN disk, or surrounded byan accretion disk can be such examples [23, 24]. In thissection, we examine this example.Now we give a demonstration of the explained methodfor three kinds of messengers; gravitational waves (GWs),high-energy neutrinos (neutrinos hereafter), and gamma-ray bursts (GRB). We assume a model with continuoussingle emissions for each messenger type, i.e. no repeatedor periodic emission for multi-messenger or single mes-senger emissions. As mentioned before, there are searchesfor multi-messenger detection of all the three combina-tions of two of these messengers [25–27]; but there is notriple messenger search. In this search, the start andend times of the GW emission or the gamma-ray emis-sion can be estimated well due to having a continuousdetection amplitude, although for neutrinos it is hardto estimate when the emission starts or ends since upto now no continuous cosmic high-energy neutrino fluxhas been detected. High-energy neutrino emissions aredetected in low numbers, generally as a single neutrino.For GWs and GRBs, the detection decision is essentiallybased on the detected continuous total energy, whereasfor neutrinos, it is based on each neutrino’s characteris-tics. Therefore it is more appropriate to separate our sig-nal sub-hypotheses based on different detected neutrinocounts (including all the characteristics), i.e., a coincidentGW–GRB– n neutrino detection. We will denote our sub-hypotheses with the notation H s = { GW,GRB,a } ,s = { b } , ..., where the sets s i in the subscript represent detectionsfrom different astrophysical sources. If the set has GW or GRB in it, that means GW or GRB emission was de-tected from that source. Finally, the sets include positiveintegers (for example, a , b ), which represent the numberof detected high-energy neutrinos from each source.For concreteness of the example we consider the groundbased interferometric detectors such as LIGO [28], Virgo[29] or KAGRA [30] for GWs, IceCube [31] for neutrinosand F ermi [32, 33] for GRBs. The detection outcomesfor GWs are x GW = { t GW , D , F} which are the detectiontime of the GW, its joint volume localization-isotropicequivalent emission energy estimation as a 4 dimensionalprobability distribution, and the estimated false alarmrate. If the joint volume localization-isotropic equivalentemission energy estimation is not explicitly provided; itcan be derived from the 3 dimensional volume localiza-tion, the detected signal energy, and the antenna pat-tern at the time of the detection. The detection out-comes for high energy neutrinos are x ν = { t ν , Ω ν , σ ν , (cid:15) ν } ,which are the detection times of the neutrinos, their ex-pected sky positions, the angular errors on the sky lo-calizations, and their reconstructed energies. The lo-calization of neutrinos are approximated as two dimen-sional Gaussian distribution, and the angular error cor-responds to one standard deviation [34]. The detectionoutcomes for GRBs are x γ = { t γ , Ω γ , κ γ , S , E} , whichare the detection time of the GRB, its expected sky po-sition, the angular error on the sky localization, its p -value, and the estimated detected energy. The localiza-tion of gamma rays is approximated as von Mises-Fisherdistribution (which is also called Gaussian distributionon the sphere). The angular error corresponds to theradius of the circle containing 68% of the probability(note that this is different than the standard deviation ofa two dimensional Gaussian distribution which contains ∼
49% of the probability) [35]. Our source parametersare θ = { r s , Ω s , t s , E GW , E ν , E γ } , which are the distanceof the source, its sky position, the retarded time of theevent, and the isotropic equivalent emission energies inGWs, high-energy neutrinos, and gamma rays. The com-plete model includes the emission delays of the messen-gers and the source rates as well, which are explainedthroughout when they are used.We will first write down the detection likelihoods inEq. (7) which encompasses the ones in Eqs. (4) and(5). For a short notation we will denote the sets of de-tection outcomes from one source { x ji , x ki , ... } as x i foreach detector. A. Detection likelihoods
We start with the GWs. The signal likelihood can beexpanded as P ( x GW | θ , H s ) = P ( t GW , D , F| t s , r s , Ω s , E GW , H s )= P ( t GW | t s , H s ) P ( F| t GW , r s , Ω s , E GW , H s ) × P ( D| t GW , F , r s , Ω s , E GW , H s ) (9) The temporal distribution of t GW is assumed to beuniform around t s : P ( t GW | t s , H s ) = ( t + GW − t − GW ) − for t GW − t s ∈ [ t − GW , t + GW ] and 0 otherwise. We take − t − GW = t + GW = 250 s as in [17, 36].By using the Bayes’ rule, we expand the likelihood forvolume localization. P ( D| t GW , F , r s , Ω s , E GW , H s )= P ( r s , Ω s , E GW | t GW , D , F , H s ) P ( D| t GW , F , H s ) P ( r s , Ω s , E GW | t GW , F , H s ) (10)The first term in the numerator is the D distributionitself. We use Bayes’ rule for the denominator to havethe form P ( D| t GW , r s , Ω s , E GW , F , H s )= P ( F| t GW , H s ) P ( r s , Ω s , E GW | t GW , H s ) P ( F| t GW , r s , Ω s , E GW , H s ) × D ( r s , Ω s , E GW ) P ( D| t GW , F , H s ) (11)The term P ( F| t GW , H s ) can be computed by integrat-ing the likelihood for fixed source parameters (which canbe obtained from calculations or simulations) over thesource parameters. P ( F| t GW , H s ) = (cid:90) P ( F| t GW , r s , Ω s , E GW , H s ) × P ( r s , Ω s , E GW | t GW , H s ) dr s d Ω s dE GW (12)Similarly, for the null hypotheses we expand the likeli-hood. P ( x GW | H n ) = P ( t GW | H n ) P ( F| t GW , H n ) × P ( D| t GW , F , H n ) (13) P ( F| t GW , H n ) can be found empirically, i.e. throughthe unphysical time shifted coincidences. We assumethe terms P ( D| t GW , F , H s ) and P ( D| t GW , F , H n ) donot depend on the hypotheses and are equal to eachother, hence cancel in the overall expression. Finallythe noise triggers are assumed to be Poisson events andhence can uniformly occur in the observation period T obs , P ( t GW | H n ) = T − obs . We note that at the end of the fullcalculation, the end result does not depend on T obs ; butwe do not drop it throughout for clarity.Next we move on the signal likelihoods for neutrinosand expand similarly. P ( t ν , (cid:15) ν , σ ν , Ω ν | θ , H s ) = P ( t ν | t s , H s ) P ( (cid:15) ν | Ω s , H s ) × P ( Ω ν , σ ν | (cid:15) ν , Ω s , H s ) (14)The temporal distribution of t ν is also assumed to beuniform around t s : P ( t ν | t s , H s ) = ( t + ν − t − ν ) − for t ν − t s ∈ [ t − ν , t + ν ] and 0 otherwise. We take − t − ν = t + ν = 250 s asin [17, 36].The estimated source localization from the detectioncan be written as P ( Ω s | (cid:15) ν , σ ν , Ω ν , H s ) = e −| Ω ν − Ω s | σ ν πσ ν (15)However we need the probability P ( σ ν , Ω ν | (cid:15) ν , Ω s , H s )which we expand with Bayes’ rule as P ( σ ν , Ω ν | (cid:15) ν , Ω s , H s )= P ( Ω s | (cid:15) ν , σ ν , Ω ν , H s ) P ( σ ν , Ω ν | (cid:15) ν , H s ) P ( Ω s | (cid:15) ν , H s )= e −| Ω ν − Ω s | σ ν πσ ν P ( σ ν , Ω ν | (cid:15) ν , H s ) P ( Ω s | (cid:15) ν , H s )= P ( σ ν , Ω ν | (cid:15) ν , H s ) e −| Ω ν − Ω s | σ ν P ( (cid:15) ν | H s )2 πσ ν P ( (cid:15) ν | Ω s , H s ) P ( Ω s | H s ) (16)By assuming a power law with exponent -2 for the energydistribution of neutrinos [37] and by using the effectivearea of the neutrino detector A eff ( (cid:15) ν , Ω s ) we write theterm P ( (cid:15) ν | H s ) as P ( (cid:15) ν | H s ) = (cid:82) A eff ( (cid:15) ν , Ω s ) (cid:15) − ν P ( Ω s | H s ) d Ω s (cid:82) (cid:15) max (cid:15) min (cid:82) A eff ( (cid:15) (cid:48) ν , Ω s ) (cid:15) (cid:48)− ν P ( Ω s | H s ) d Ω s d(cid:15) (cid:48) ν (17) (cid:15) min , (cid:15) max are 100 GeV and 100 PeV for IceCube. The P ( (cid:15) ν | Ω s , H s ) terms in Eqs. (14) and (16) cancel.Next, we expand the null hypothesis likelihood simi-larly. P ( t ν , (cid:15) ν , σ ν , Ω ν | H n ) = P ( t ν | H n ) P ( (cid:15) ν , σ ν , Ω ν | t ν , H n )= P ( t ν | H n ) P ( (cid:15) ν | t ν , H n ) P ( σ ν , Ω ν | (cid:15) ν , t ν , H n ) (18) P ( t ν | H n ) = T − obs and P ( (cid:15) ν | t ν , H n ) can be found empiri-cally from detector characteristics and past observations.The time dependency of the last term comes from the an-nual modulation due to Earth’s motion around the Sunand can be expressed with a function T ( t ν , (cid:15) ν , Ω ν ) whoseaverage over one year for every ( (cid:15) ν , Ω ν ) pair is one. P ( σ ν , Ω ν | (cid:15) ν , t ν , H n ) = P ( σ ν , Ω ν | (cid:15) ν , H n ) T ( t ν , (cid:15) ν , Ω ν )(19)The terms P ( σ ν , Ω ν | (cid:15) ν , H s ) and P ( σ ν , Ω ν | (cid:15) ν , H n ) do notdepend on the hypotheses and cancel in the overall ex-pression.Third, we move on the likelihoods for GRBs and ex-pand similarly. P ( t γ , Ω γ , κ γ , S , E| θ , H s ) = P ( t γ | t s , H s ) × P ( S| t γ , Ω s , r s , E γ , H s ) P ( Ω γ , κ γ , E|S , Ω s , t γ , E γ , H s )(20)The temporal distribution of t γ is also assumed to beuniform around t s : P ( t γ | t s , H s ) = ( t + γ − t − γ ) − for t γ − t s ∈ [ t − γ , t + γ ] and 0 otherwise. We take t − γ = 100 s and t + γ = 250 s from [36]. The source localization anddetected energy from the detection can be summarisedas P ( Ω s , r s , E γ | t γ , Ω γ , κ γ , E , S , H s )= κ γ πsinh ( κ γ ) e κ γ cos ( | Ω γ − Ω s | ) δ ( E − η E γ πr s ) (21)where η is a constant describing the detection ef-ficiency of the detector. Similarly to the neutrinocase we need a different but a related probability P ( Ω γ , κ γ , E|S , Ω s , r s , t γ , E γ , H s ) which we expand viathe Bayes’ rule. P ( Ω γ , κ γ , E|S , Ω s , r s , t γ , E γ , H s )= P ( Ω s , r s , E γ |S , Ω γ , κ γ , E , t γ , H s ) × P ( Ω γ , κ γ , E|S , t γ , H s ) P ( Ω s , r s , E γ |S , t γ , H s )= P ( Ω s , r s , E γ |S , Ω γ , κ γ , E , t γ , H s ) P ( S| Ω s , r s , E γ , t γ , H s ) P ( Ω s , r s , E γ | t γ , H s ) × P ( Ω γ , κ γ , E|S , t γ , H s ) P ( S| t γ , H s ) (22)The P ( S| Ω s , r s , E γ , t γ , H s ) terms in Eqs. (20) and (22)cancel. The term P ( S| t γ , H s ) can be computed by inte-grating the likelihood for fixed source parameters (whichcan be obtained from calculations or simulations) overthe source parameters. P ( S| t γ , H s ) = (cid:90) P ( S| Ω s , r s , E γ , t γ , H s ) × P ( Ω s , r s , E γ | t γ , H s ) d Ω s dr s dE γ (23)We expand the null hypothesis likelihoods as P ( t γ , Ω γ , κ γ , S , E| θ , H n ) = P ( t γ | H n ) × P ( S| t γ , H n ) P ( Ω γ , κ γ , E|S , t γ , H n ) (24) P ( Ω γ , κ γ , E|S , t γ , H s ) and P ( Ω γ , κ γ , E|S , t γ , H n ) termsdo not depend on hypotheses and cancel in the overallexpression. P ( S| t γ , H n ) is equal to the p value of theevent and P ( t γ | H n ) = T − obs . B. Prior subhypothesis probabilities
Now we move on the prior probabilities for each sub-hypothesis. These are found by assuming each detectioncandidate trigger (noise or astrophysical origin) is a Pois-son event. The expected counts for the Poisson processesare found by the known noise trigger rates R bg,ξ and the assumed true astrophysical source rates ˙ n trueξ for the mes-senger ξ . We are interested in the observable source rates˙ n ξ for GWs and GRBs which have detection cuts in termsof the signal to noise power ratio or photon count. We de-fine ρ ( E GW r s , Ω s , t s ) and I ( E γ r s , Ω s , t s ) functions as the cutfunctions and the detection thresholds ρ th and I th . ρ canbe taken as the network signal-to-noise ratio for GWs and I as the detected intensity or the photon count. Thosefunctions take in account the effective antenna patternof the GW detector network (by accounting the differ-ent sensitivities of the detectors too) and the view of the F ermi satellite. Furthermore, we assume beaming forneutrino and gamma-ray emission from the same open-ing with a beaming factor f b ( ∼ − n GW = (cid:90) ˙ n trueGW P ( θ )[ ρ ( E GW r s , Ω s , t s ) ≥ ρ th ] d θ (25)The binary bracket notation [ ζ ] is 1 if ζ is true and 0 iffalse. For a GRB only source˙ n γ = f − b (cid:90) ˙ n trueγ P ( θ )[ I ( E γ r s , Ω s , t s ) ≥ I th ] d θ (26)For a multi-messenger source it is˙ n GW,γ = f − b (cid:90) ˙ n trueGW,γ P ( θ )[ ρ ( E GW r s , Ω s , t s ) ≥ ρ th ] × [ I ( E γ r s , Ω s , t s ) ≥ I th ] d θ (27)For neutrinos we are interested in the observable neutrinorate rather than the observable source rate.˙ n ν = f − b (cid:90) ˙ n trueν P ( θ ) (cid:104) n ν ( E ν , r s , Ω s ) (cid:105) d θ (28) (cid:104) n ν ( E ν , r s , Ω s ) (cid:105) is the detector specific expected numberof detected neutrinos from a source with given locationand emission energy which scales linearly with E ν r s anddepends on the effective area. For a multi-messenger de-tection with neutrinos the interesting quantity would be˙ n GW,ν,γ = f − b (cid:90) ˙ n trueGW,ν,γ P ( θ ) (cid:104) n ν ( E ν , r s , Ω s ) (cid:105)× [ ρ ( E GW r s , Ω s , t s ) ≥ ρ th ][ I ( E γ r s , Ω s , t s ) ≥ I th ] d θ (29)For clarity, let’s demonstrate a specific subhypothesis H s = { GW,GRB,a } ,s = { b } for detected α GWs, β GRBs and µ neutrinos in total. As a reminder, that subhypoth-esis corresponds to signal detections from two sources;from the first one, s , a GW, a GRB and a neutrinosare detected, and from the second one, s , only b neu-trinos are detected. We will denote the probability ofoccurrence of d Poisson events with an expectation λ as P oi ( d, λ ) = λ d e − λ d ! . In this subhypothesis the first sourceis clearly a multi-messenger source; but the second onecan be a multi-messenger source from which the GW orthe GRB or both were not detected, or it can be simplya source which only emits neutrinos. We consider all ofthese 4 possible cases. P ( H s = { GW,GRB,a } ,s = { b } | θ , θ ) = P oi ( µ − a − b, R bg,ν T obs ) P oi ( α − , R bg,GW T obs ) P oi ( β − , R bg,γ T obs ) × P oi ( a, (cid:104) n ν ( E ν , r s , Ω s ) (cid:105) )[ ρ ( E GW r s , Ω s , t s ) ≥ ρ th ][ I ( E γ r s , Ω s , t s ) ≥ I th ] P oi ( b, (cid:104) n ν ( E ν , r s , Ω s ) (cid:105) ) × { P oi (2 , ˙ n GW,ν,γ T obs ) P oi (0 , ( ˙ n GW − ˙ n GW,ν,γ ) T obs ) P oi (0 , ( ˙ n ν − ˙ n GW,ν,γ ) T obs ) P oi (0 , ( ˙ n γ − ˙ n GW,ν,γ ) T obs ) × [ ρ ( E GW r s , Ω s , t s ) < ρ th ][ I ( E γ r s , Ω s , t s ) < I th ]+ P oi (1 , ˙ n GW,ν,γ T obs ) P oi (1 , ˙ n GW,ν T obs ) P oi (0 , ( ˙ n GW − ˙ n GW,ν − ˙ n GW,ν,γ ) T obs ) × P oi (0 , ( ˙ n ν − ˙ n GW,ν − ˙ n GW,ν,γ ) T obs ) P oi (0 , ( ˙ n γ − ˙ n GW,ν,γ ) T obs )[ ρ ( E GW r s , Ω s , t s ) < ρ th ]+ P oi (1 , ˙ n GW,ν,γ T obs ) P oi (1 , ˙ n ν,γ T obs ) P oi (0 , ( ˙ n GW − ˙ n GW,ν,γ ) T obs ) P oi (0 , ( ˙ n ν − ˙ n ν,γ − ˙ n GW,ν,γ ) T obs ) × P oi (0 , ( ˙ n γ − ˙ n ν,γ − ˙ n GW,ν,γ ) T obs )[ I ( E γ r s , Ω s , t s ) < I th ]+ P oi (1 , ˙ n GW,ν,γ T obs ) P oi (0 , ( ˙ n GW − ˙ n GW,ν,γ ) T obs ) P oi (1 , ( ˙ n ν − ˙ n GW,ν,γ ) T obs ) P oi (0 , ( ˙ n γ − ˙ n GW,ν,γ ) T obs ) } (30) C. Source parameter distributions
Finally, we explain the required distributions for sourceparameters. First, we write the complete distribution.The sources are distributed such that the event rate is uniform in the comoving spacetime. There can be manydifferent models for the emission energies. Here we pro-vide only a naive example. We assume log uniform dis-tributions for GW, neutrino, and GRB emission energies[17]. We take the limits of neutrino and GRB emissionsto be 10 − ergs [26, 38, 39] and GW limits to bebetween 0 . − M (cid:12) c (assuming ∼
5% of the mass isemitted in a merger [40]). The event time is distributeduniformly in the observation time. P ( θ ) = r s π (1 + z ( r s )) E GW E ν E γ T obs N r log (100) (31) N r is the normalization constant for r s . z ( r s ) is the red-shift and the factor (1 + z ( r s )) in the denominator ac-counts for the dilution of sources in space and the timedilation due to Hubble expansion.There are three conditional source distributions usedin the likelihoods. The one in the GW part is P ( r s , Ω s , E GW | t GW , H s ) = r s π [ ρ ( E GW r s , Ω s , t GW ) ≥ ρ th ](1 + z ( r s )) E GW N (cid:48) r log (100)(32) N (cid:48) r is the normalization constant. Here we ignored theeffect of using t GW instead of t s in the ρ function. Forgreater accuracy a new function can also be defined.The conditional distribution in the neutrino part is P ( Ω s | H s ) = (cid:82) (cid:15) max (cid:15) min A eff ( (cid:15) ν , Ω s ) (cid:15) − ν d(cid:15) ν (cid:82) (cid:82) (cid:15) max (cid:15) min A eff ( (cid:15) ν , Ω (cid:48) s ) (cid:15) − ν d(cid:15) ν d Ω (cid:48) s (33)The conditional distribution in the GRB part is P ( r s , Ω s , E γ | t γ , H s ) = r s [ I ( f b E γ r s , Ω s , t γ ) ≥ I th ]4 π (1 + z ( r s )) E γ N (cid:48)(cid:48) r log (100)(34) N (cid:48)(cid:48) r is the normalization constant. Here we also ignoredthe effect of using t γ instead of t s in the I function. Forgreater accuracy a new function can also be defined.With the guidance provided in this section, a realtimemulti-messenger search for GWs, neutrinos and GRBscan be constructed. V. CONCLUSION
In this paper, we addressed the problem of optimalmulti-messenger searches. Having more messengers willnot only make us better understand their sources; butalso can increase the significance of sub-threshold singlemessenger detections and increase the rate of detectionswithout a necessary upgrade to the detectors.We showed that a model-independent optimal solutiondoes not exist. We provided a Bayesian solution thatis scalable to any number of messengers. It is based onconstructing a test statistic by combining different sub-hypotheses via using their predicted rates according toa model. This gives the highest power for the regularfrequentist hypothesis test for the assumed model. As aBayesian solution, this method’s performance is depen-dent on the accuracy of the current models. The de-scribed method is completely scalable and applicable toany number and any kind of messengers.Finally, we examined the use case for a search forjoint GW-neutrino-GRB emissions. Although there aresearches for all the three combinations of two of thesemessengers [25–27], this is the first examination of thetriple messenger search, which can be applied in real timee.g., similarly to [41, 42].
ACKNOWLEDGMENTS
The authors thank Benjamin Farr and Gregory Ashtonfor useful feedback. The authors are thankful for the sup-port of Columbia University in the City of New York andthe National Science Foundation grant PHY-2012035.D.V. acknowledges Jacob Shaham Fellowship. I.B. ac-knowledges support from the National Science Founda-tion under grant PHY-1911796 and the Alfred P. SloanResearch Foundation. This document was reviewed bythe LIGO Scientific Collaboration under the documentnumber P2000377. [1] G. Magli, Nexus Netw J , 337–346 (2016).[2] M. Hoskin,