A Method for an Untriggered, Time-Dependent, Source-Stacking Search for Neutrino Flares
AA Method for an Untriggered, Time-Dependent,Source-Stacking Search for Neutrino Flares
The IceCube Collaboration † The IceCube Collaboration ∗ http://icecube.wisc.edu/collaboration/authors/icrc19_icecube E-mail: [email protected], [email protected] Recent results from IceCube regarding TXS 0506+056 suggest that it may be useful to test thehypothesis of multiple neutrino flares, where each flare is not necessarily accompanied by acorresponding gamma-ray flare. An untriggered, time-dependent, source-stacking search wouldbe optimal for testing this hypothesis, however such an analysis has yet to be applied to the fullduration of IceCube data. Here, we discuss one possible way of constructing such an analysis,with a global test statistic obtained by summing the individual test statistics associated witheach signal-like flare in the sample. This has the additional advantage of assessing the timestructure of each source when running the analysis over a source catalog. We show that, fora signal consisting of many small flares, this style of analysis represents a significant increasein discovery potential over both the existing single-flare fit and the time-integrated stackingmethods. Potential source catalogs are examined in combination with this method, including thepossibility of a “self-triggered” catalog consisting of the locations of the highest energy northernsky events in the IceCube sample.
Corresponding authors:
William Luszczak † , Jim Braun , Albrecht Karle Dept. of Physics and Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin,Madison, WI 53706, USA36th International Cosmic Ray Conference -ICRC2019-July 24th - August 1st, 2019Madison, WI, U.S.A. ∗ For collaboration list, see PoS(ICRC2019) 1177. † Speaker. c (cid:13) Copyright owned by the author(s) under the terms of the Creative CommonsAttribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND 4.0). http://pos.sissa.it/ a r X i v : . [ a s t r o - ph . H E ] A ug ultiflare Stacking William Luszczak
1. Introduction
Recent results from the IceCube Neutrino Observatory suggest the possibility of multiple un-triggered neutrino flares from the same object [1][2]. The method used in [1] tests the hypothesisof multiple flares for TXS 0506+056 by combining the single, most significant flares found overdifferent pre-defined data taking periods. While this method works well in the case that the data isdominated by a single, significant neutrino flare, it is not sensitive to the case of multiple, similarly-sized flares produced by the same source. Here, we present a method specifically tailored to testingthe hypothesis that a particular source flares many (>1) times. We will refer to this new method as"multiflare stacking". The primary difference between multiflare stacking and the existing methods[3], is that multiflare stacking attempts to fit all flares in the sample, not merely the largest flare.This method is particularly well suited for testing the hypothesis of many, moderate-to-weaklyflaring sources.
2. Description of the method
Conceptually, multiflare stacking can be broken down into the following steps:1. Create a set of "test windows", defined by a start and stop time. This is done by taking allpossible pairs of events that pass a pre-defined weight threshold.2. Calculating a test statistic associated with each test window, describing the hypothesis thatthat particular window was a flare produced by the source. Discard all windows with teststatistic below some threshold: These windows are "background-like", while the remainingwindows (with test statistic > this threshold) are "signal-like".3. For the remaining signal-like windows, remove windows that overlap with a window withhigher test statistic. Sum the remaining
T S j to obtain a multiflare test statistic, (cid:102) T S , associatedwith a particular source.4. (For the case of considering a source catalog) Do a binomial test to optimize for the numberof sources to stack together.
For the purpose of multiflare stacking, we will restrict ourselves to only considering flareswhich begin and end on an event. This means that for a data sample containing N events, thereare at most N ( N − ) / N ( N − ) / S / B ratio greater than some threshold will be considered for forming test windows. Here, S / B refers to the ratio of the spatial and energy components of the signal and background PDFsdefined in equations 4 and 8 of [3]. This significantly reduces the number of test windows to amore manageable number. This process is shown in Figure 1. This S / B event selection is onlyused for the test window determination. For the actual analysis all N events will be used.2 ultiflare Stacking William Luszczak
Figure 1: The process of selecting test windows for which a test statistic will be calculated. Testwindows may overlap, and windows are only formed from events passing the S / B threshold. Notethat while the S / B threshold is applied to determine which windows will be considered, all events(regardless of S / B value) will be used when calculating a test statistic. For each of the test windows defined in section 2.1, we calculate a flare test statistic (
T S j associated with that specific flare hypothesis (that the flare starts at the time of the first seed event,and ends at the time of the second seed event). We identify individual test windows with theindex j . This flare test statistic is of the form presented in section 3 of [3], with the modificationthat we use a box-shaped flare hypothesis instead of a gaussian. . Signal PDFs ( S i j where theindex i indicates the event) are defined as in equations 2.1-2.3. Events (identified by the index i )are defined by a arrival direction ( (cid:126) r i ), energy ( E i ), and time ( t i ). The signal PDF is composed ofa spatial term ( R i j ( (cid:126) r i ) ), defined by the location (cid:126) r i and event angular error σ i . The energy PDF, E ( E i j | γ j ) , describes the probability of obtaining an event energy ( E i ) from a signal event producedby a flare with spectrum γ j , and is a function both of γ j and detector effective area. Note that weinclude a j index with the flare spectrum, γ j , as this method will fit the spectrum of each flareindividually. In the definition of the temporal PDF (equation 2.3), H ( t ) , is the Heaviside function, ∆ T j is the test window duration, and t startj and t endj are the start and end time of the j th test window. S i j = R i j ( (cid:126) r i ) × E ( E i | γ j ) × T i j ( t i ) (2.1) R i j ( (cid:126) r i ) = πσ i exp [ − r i σ i ] (2.2) T i j ( t i ) = H ( t endj − t i ) H ( t i − t startj ) ∆ T j (2.3)The background PDF is defined in equation 2.4. Similar to the signal PDFs, the backgroundenergy and spatial PDFs are functions of the detector effective area, while the temporal backgroundPDF, ∆ T data , is simply a uniform distribution over the entire data taking period. B i = R bg ( (cid:126) r i ) E ( E i | Atm ν ) ∆ T data (2.4)3 ultiflare Stacking William Luszczak
We can then define our test window likelihood and associated test statistic, shown in equations2.5 and 2.6. Here, N refers to the total number of events in the sample (without the S / B requirementmentioned above), and n s j is the number of signal events contained within the j th flare. The flarelikelihood (equation 2.5) can be minimized with respect to n s j and γ j to obtain the flare test statisticin equation 2.6, where ˆ n s j and ˆ γ j refer to the best-fit values resulting from minimizing equation2.5. This test statistic describes the hypothesis that the j th flare is in fact a signal flare. Unlikethe standard method, we do not minimize our likelihood with respect to the flare parameters ∆ T j , t endj or t startj , as instead we seek to describe the hypothesis of a flare with these parameters fixed tospecific values. L j ( n s j , γ j ) = N ∏ i = ( n s j N S i j + ( − n s j N ) B i ) (2.5) T S j = − (cid:20) ∆ T data ∆ T j × L j ( n s j = )) L j ( n s j = ˆ n s j , γ j = ˆ γ j ) (cid:21) (2.6)The ∆ T data / ∆ T j term is intended to correct for the fact that there are many small windows thatformed in section 2.1, so we penalize their test statistic to account for the additional trials associatedwith short test windows. Note that while test windows are defined by events that pass the thresholdin step 1, the i index in equation 2.5 runs over all events, not only those which pass the S / B cut.This process can then be repeated for every test window defined in section 2.1. The end resultof this is a list of flares (start/stop times), and corresponding T S j ’s. We then discard all windowswith T S j <
0. These windows are "background-like", while the remaining windows, with
T S j > Individual events should not be able to contribute to more than one flare. To account for this,we employ the following decorrelation procedure: test windows are ordered by descending
T S j (calculated above). Iterating down the list, discard all test windows that overlap in time with adifferent test window with a larger T S j . The result of this step is a list of completely decorrelated(non-overlapping) test windows. This list can be interpreted as a "neutrino flare curve", describingthe set of neutrino flares associated with a particular source location. Figure 2 shows an exampleflare curve generated by this procedure.We can now define the "global" test statistic, (cid:102) T S to be the sum of all the
T S j remaining afterthis decorrelation procedure, as in equation 2.7, where the sum in equation 2.7 is only over non-overlapping flares. The best fit value for the number of signal events, ˜ n s , can be obtained in asimilar manner by simply summing up the best fit ˆ n s associated with all the non-overlapping flares. (cid:102) T S = ∑ j , T S j > T S j (cid:101) n s = ∑ j , T S j > ˆ n s j (2.7)Figure 3 shows the best fit ˜ n s versus the injected value for both this method and the traditionalsingle flare fit. For a signal composed of many flares, multiflare stacking recovers the injectedsignal much more accurately. 4 ultiflare Stacking William Luszczak
Figure 2: Multiflare stacking applied to a single source (a) Top: Events that pass the S / B threshold cutfor generating windows, mentioned in section 2.1,for a source located at (ra, dec) = (77 . ◦ , 5 . ◦ ).This source contains three injected flares, centered atMJD=55246.3, 55807.4, and 57632.7. Bottom: Theresultant flare curve calculated by applying multiflarestacking to this source. The single flare method fitsonly the flare at MJD = 55246.3 (with TS = 22.66),while the multiflare method fits all flares together,with a global test statistic of (cid:102) T S = 60.52. (b) The background test statistic distributions for asingle source, located at declination = 5 . ◦ . Thebackground multiflare test statistic distribution isshown in blue, while the single flare method is shownin orange. The vertical lines represent the test statis-tics associated with the flare curve in figure 2a. Thesingle flare p-value for this source is 4 . × − (3.29 σ ), while the multiflare p-value is 8 . × − (4 . σ ). Figure 3: Fit versus injected values for the number of signal events ( n s ), for both the existing singleflare method (left) and the multiflare method presented in these proceedings (right). In both casesa signal composed of many flares was injected, with flare intensities drawn according to the powerlaw distribution introduced in section 3. Unsurprisingly, multiflare stacking is much more capableof recovering the number of injected signal events for a signal of this form. We can easily extend this process to a source catalog, rather than a single source. The global5 ultiflare Stacking
William Luszczak test statistic is simply given by equation 2.8, where here the sum runs over all sources in the sourcecatalog. (cid:102)
T S all = ∑ sources (cid:102) T S (2.8)For large source catalogs, if only a few sources are "truly" flaring, they may be overwhelmedby the large number of background sources. This can be remedied by "binomial stacking" ofsources. Instead of stacking the multiflare test statistic for all sources (equation 2.8), we insteademploy the following procedure to optimize the number of sources to stack ( k ). For a given dataset, the sources are ordered by their multiflare test statistic, (cid:102) T S . A p-value for k = , , , ... N srcs iscalculated, and subsequently the k that produces the minimum p-value is selected. An additionaltrial factor is attached from a newly created background test statistic distribution, correspondingto the number of stacked sources, similar to the process described in [4]. The additional powerobtained by stacking in this manner outweighs the penalty associated with the trial factor, resultingin an increased discovery potential (see section 4). (cid:102) T S all = k best ∑ m = (cid:102) T S m (2.9)
3. Characterizing an Injected Signal
In order to test the performance of multiflare stacking, we require a description of the intensitydistribution of observed flares (i.e. how many events are in each flare, and how many flares thereare total). For flares of intensity I , the number of flares with that intensity N ( I ) can be described bya power law with a lower cutoff (equation 3.1). This is primarily motivated by the distribution ofgamma ray flares found in [5]. Note, however, that this is a descriptive power law, not a predictive power law. The parameters here are intended as a flexible way to describe a wide array of whatthe distribution of flares may actually look like. We make no claims as to what these parametersactually are beyond the existing limits calculated in section 4. N ( I ) = A m ( α − I o )( II o ) − α (3.1)This power law is specified by three parameters: • A m describes the overall power law normalization (how many flares exist in the sample) • I o describes the lower intensity cutoff. This is a parameter that scales with the flux, butroughly corresponds to the minimum flare size. • α describes the number of larger versus smaller flares. A larger value of α corresponds toproportionally more smaller flares, while a smaller α will produce more larger flares.
4. Catalogs and Discovery Potential
In this section we will explore the 3 σ discovery potential of multiflare stacking applied totwo different catalogs. We define the 3 σ discovery potential to be the flare intensity distributionnormalization ( A m ) required for 50% of the injected signal multiflare test statistics to be higherthan the 3 σ threshold calculated from the background multiflare test statistic distribution. Under6 ultiflare Stacking William Luszczak this construction, the 3 σ discovery potential is a surface in the space spanned by possible valuesof A m , α and I o . For plots in this section, we have chosen to show slices of this space with fixed α ( α = A m as a function of I o . Since we have yet to observea flaring neutrino source, the allowed values of α are unconstrained, we have simply chosen aspecific value to illustrate a region where multiflare stacking is relevant. This is purely for thepurposes of demonstration, as this method has significant discovery potential for a wide range ofinjected values of α . Due to previous results of time-integrated stacking searches applied to Fermi 2LAC blazars[6], it may be beneficial to test the hypothesis of multiple flares for the catalog of Fermi 3LACblazars. We select the top 500 flux sources in this catalog purely due to computational limitations,but otherwise include no additional cuts.The currently existing limits for this catalog are plotted in blue in figure 4a, and can be cal-culated analytically by integrating equation 3.1. The "single flare limits" correspond to performinga binomial test on the Fermi 3LAC blazar catalog (a 3 sigma discovery corresponds to roughly 5pre-trial 3-sigma injections), as in [4]. The "time integrated limits" correspond to the conditionthat there are not more events in the sample than are allowed by the existing time-integrated fluxlimits reported in [6]. (a) The 3 σ multiflare discovery potential for a cata-log of Fermi 3LAC blazars. The dashed line corre-sponds to the limits set in [6]. (b) The existing limits and 3 σ discovery potential forthe "self-triggered" catalog described in section 4.2.The dashed line in this plot refers to the diffuse fitdone in [7]. Figure 4: Multiflare discovery potential. Injected flare durations are fixed to 100 days, and theinjected events follow an E − . spectrum. Note that the discovery potential improves a small, butsignificant amount with the inclusion of binomial-style stacking (section 2.4). In both figures, thesolid blue line corresponds to the observation of a single, 100 day flare with 13 events.Figure 4a shows the 3 σ discovery potential for applying multiflare stacking to this source cat-alog. Note that the binomial-style stacking method proposed in section 2.4 significantly improvesthe discovery potential, improving it beyond the currently existing limits set by [6]. Also note that7 ultiflare Stacking William Luszczak for large values of I o , the discovery potential of the single flare fits is expected to exceed that ofmultiflare stacking, since this region corresponds to having only 1-2 large flares in the sample. Forvery small values of I o , the time integrated method becomes superior, since in this region "flares"are often merely a single event. Multiflare stacking is primarily sensitive to the middle region,where there are a moderate amount of medium-sized flares. Based on the diffuse fits to the atmospheric and astrophysical spectrum in [7], we can developa "self-triggered" catalog of source locations consisting of the highest energy events in the sample.This is similar in spirit to the catalog considered in [8]. For this catalog, we test the locations ofall events in the sample with energy proxy > 200 TeV. This threshold is chosen because above 200TeV, we expect 50% of events to be astrophysical in origin [7]. These "source" events are removedfrom the sample prior to calculating a test statistic to avoid signal contamination.This catalog consists of 32 high-energy events, though notably does not include IC170922A,as this data sample only extends up to 2016. Events are well localized, with the majority of eventshaving reconstructed angular error < 1 ◦ , so these locations can be considered to be point-like sourcecandidates.The discovery potential for this catalog is shown in figure 4b. The benefit of the binomialtest is smaller in this case, due to the size of the source catalog, however multiflare stacking stillrepresents a significant improvement over both the single flare and time-integrated methods.
5. Concluding Remarks
We have presented a method for testing the hypothesis that neutrino sources flare multipletimes ("multiflare stacking"). This is done by combining information from fits to all flares in thesample. Using an injected signal characterized by a power law, we explored the discovery potentialfor this method for two source catalogs: Fermi 3LAC blazars and a "self-triggered" catalog com-posed of the locations of high energy IceCube events observed between 2009 and 2016. For bothof these catalogs, multiflare stacking represents a significant improvement over existing methodsfor a signal composed of many moderately-sized flares.
References [1]
IceCube
Collaboration, M. G. Aartsen et al.,
Science (2018) 147–151.[2]
IceCube
Collaboration, M. G. Aartsen et al.,
Science (2018) eaat1378.[3] J. Braun, M. Baker, J. Dumm, C. Finley, A. Karle, T. Montaruli,
Astroparticle Physics (2010)175–181.[4] IceCube
Collaboration, E. O’Sullivan,
PoS(ICRC2019)973 (these proceedings).[5] M. Giomi, et al.,
The Astrophysical Journal Supplement (2016).[6]
IceCube
Collaboration, M. G. Aartsen et al.,
ApJ (2017) 45.[7]
IceCube
Collaboration, M. G. Aartsen et al.,
ApJ (2016) 3.[8]
IceCube
Collaboration, M. Karl,
PoS(ICRC2019)929 (these proceedings).(these proceedings).