Methods for point source analysis in high energy neutrino telescopes ∗ Jim Brauna, , Jon Dumma, Francesco De Palmaa, Chad Finleya, Albrecht Karlea, Teresa Montarulia,b 8 aUniversity of Wisconsin, Chamberlin Hall, Madison, Wisconsin 53706, USA 0 0 bOn leave of absence from Universit`a di Bari and INFN, 70126 Bari, Italy 2 n a J Abstract 0 1 Neutrino telescopes are moving steadily toward the goal of detecting astrophysi- ] cal neutrinos from the most powerful galactic and extragalactic sources. Here we h describe analysis methods to search for high energy point-like neutrino sources us- p - ing detectors deep in the ice or sea. We simulate an ideal cubic kilometer detector o based on real world performance of existing detectors such as AMANDA, IceCube, r t and ANTARES. An unbinned likelihood ratio method is applied, making use of the s a point spread function and energy distribution of simulated neutrino signal events [ to separate them from the background of atmospheric neutrinos produced by cos- 1 mic ray showers. The unbinned point source analyses are shown to perform better v than binned searches and, depending on the source spectral index, the use of energy 4 0 information is shown to improve discovery potential by almost a factor of two. 6 1 . 1 0 8 0 1 Introduction : v i X r With the construction of IceCube at the South Pole [1] and of ANTARES in a the Mediterranean Sea [2], together with existing R&D programs for a cubic kilometer array at these latitudes [3,4], neutrino astronomy is entering a very promising era. IceCube, when complete in 2011, will consist of up to 80 strings on a hexagonal grid with 124 m spacing, with each string holding 60 optical modules vertically spaced by 17 m. The instrumented part of the strings is deployed between 1450 and 2450 m deep in the ice. The experiment profits from experience acquired with the AMANDA detector, taking data since 1996 and completed in 2000 with 19 strings containing 677 total optical modules ∗ Corresponding author. Email address: [email protected] (Jim Braun). Preprint submitted to Elsevier 4 February 2008 between 1500 and 2000 m below the ice surface. Studies on the optimal config- uration for a cubic kilometer array in the Mediterranean are ongoing, and the amount of photomultipliers (PMTs) generally considered is somewhat larger than for IceCube, around 6000 [5] and up to about 9000 [6]. The community is refining methods to detect low statistics signals amongst large backgrounds. The expected background from atmospheric neutrinos in a cubic kilometer detector is of the order of 50 000 upgoing events per year after selection criteria guaranteeing good angular resolution and rejection of misreconstructed cosmic ray muons. Point-like signals of few events need to be singled out among this large number of background events. Two features distinguish signal from the background: • The angular distribution. The signal would cluster around the direction of the neutrino source (assumed here to be point-like) with a spread depending on the detector angular resolution. Angular resolution is limited by detector geometry and by the propagation characteristics of light in the medium, specifically by photon scattering and absorption. The pointing accuracy for astrophysical sources is also limited by the kinematic angle between the parent neutrino and the muon. • The energy distribution. The differential energy spectrum of the signal ex- pectedfromFermiaccelerationmechanismsisclosetoE−2,harderthanthat of atmospheric neutrinos due to the showering process in the atmosphere. The differential spectrum of atmospheric neutrinos approximately follows a power law of E−3.7 above 100 GeV. Other signatures may be used, including time dependencies of emissions mea- sured in other detectors such as correlations with gamma ray bursts or TeV gamma ray flares. However, we focus on steady emissions of neutrinos with time. The methods we have implemented exploit the two features listed above. We show that unbinned methods based on the likelihood ratio hypothesis test perform better than methods based on angular bins, and we show that the introduction of energy dependent information, e.g. the number of hit PMTs, helps discriminate signal and allows an energy spectrum reconstruction even when few events are detected on top of the background. Other unbinned meth- ods have been studied and developed by other authors [7,8,9,10]. Sec. 2 describes the simulation we use to generate realistic samples of signal and background events in a cubic kilometer detector. Sec. 3 describes the unbinned method based on a likelihood ratio analysis, comparing a signal plus background hypothesis to a background only hypothesis. This method has been applied to IceCube data for 9 strings and to 2005-6 data of the AMANDA-II detector [11,12]. In Sec. 4 we describe the performance of the methods and show results in terms of discovery potential. We also emphasize the importance of using energy related information to increase the discovery 2 16.7m 125m 1km Fig. 1. Simulated detector consisting of 81 strings, each containing 60 optical mod- ules. potential and show the ability to determine the spectral index of the neutrino source. The method is then compared to a search using angular bins, more traditionally applied in neutrino astronomy [13,14,15]. 2 Simulation of a data sample of atmospheric neutrino background and point source signal We wish to compare several neutrino point source search methods and draw general conclusions on the discovery potential for cubic kilometer scale neu- trinotelescopesunderconstructionintheiceandunderstudyintheseawater. We have performed a realistic simulation of both atmospheric neutrino events and signal events from an astrophysical neutrino point-like source. This is accomplished by a detailed detector volume simulation with simplifications expected to have negligible impact on the comparison. Thesimulateddetectorconsistsof4860opticalmodulesarrangedin81strings. The strings are evenly distributed on a square 9 × 9 grid with 125 meters separating nearest neighbors, and the modules on each string are vertically separated by 16.67 m, shown in Fig. 1. The modules are simulated as contain- ing a downward looking 10 inch photomultiplier with 20% quantum efficiency. A multiplicity trigger requires that at minimum 14 modules register a pho- 3 ton from an event. We neglect trigger time windows and photon hits from photomultiplier dark noise, since hits unrelated to the track can presumably be removed with coincidence requirements, filtering strategies, and topological cuts. We assume the detector is located at the South Pole. A different location would imply a different visibility of sources in the sky, since at the South Pole half of the sky is always visible while the other half is inaccessible. At other latitudes the detector is always blind to less than one half of the sky, always sensitive to a region of similar size in the opposite hemisphere, and sensitive to the remainder of the sky for a fraction of the day. We have prepared an algorithm to simulate the muon neutrino interactions in a volume and propagate the secondary muon in the ice, producing the light detected by the PMTs. For neutrino generation, we use an updated versionofthesimulationdescribedin[16].WeusethemorerecentCTEQ6[17] structure functions for the deep inelastic cross section of muon neutrinos and simulate the Earth density profile in [18] to account for the absorption of high energy neutrinos. The muon is propagated to the instrumented region using the MUM propagation code [19]. We model the muon energy loss within the detector as dE/dx = a + bE, where a = 0.268 GeV/m accounts for ionization energy losses, and b = 0.00047 m−1 for bremsstrahlung, pair production, and photonuclear interactions [20]. The number of PMTs recording light, or ‘hit’, increases with the energy loss rate, and thus the energy, of the muon. Since the energy dependent processes are stochastic, treating them as continuous may overestimate the number of modules hit by a small amount at PeV energies. However, high energy muons where this effect is significant are very bright and tend to produce photon hits in many more modules than the trigger threshold of 14, so we conclude that the impact on event triggering is negligible. Photon propagation in the detector medium has a significant impact on the detector response. We assume a homogeneous detector medium in which the photondensityatthesensorcanbedescribedasasimplefunctionwithrespect to the distance of the muon track. Photons are propagated with an effective scattering length of 21 m and an absorption length of 120 m, which are typical values for South Pole ice [21]. Fig. 2 shows the assumed photon density for a minimum ionizing muon as a function of distance for the case of ice. For a muon with arbitrary energy loss dE/dx, the photon density is scaled by the equivalent energy loss of a minimum ionizing muon, i.e. scaled by dE/dx / (0.268 GeV/m). The number of photoelectrons recorded from each PMT is a Poisson random variable with mean equal to the product of this photon den- sity, PMT photocathode area, and PMT quantum efficiency. For each PMT, a hit is determined by randomly sampling this Poisson probability of observing at least one photoelectron. We have generated 1011 upgoing neutrino events isotropically entering the Earth. The events are drawn from an E−1.4 energy spectrum between 10 and 4 103 2]m ns/102 o t o h P 10 [ y t si 1 n e D on 10-1 t o h P10-2 0 20 40 60 80 100 120 140 160 180 200 220 240 Radius from Track [m] Fig. 2. Photon density as a function of the radius from a minimum ionizing muon track in ice with an effective scattering length of 21 meters and absorption length of 120 meters. 109 GeV. This hard spectrum allows efficient generation of high energy events. Eventspassingthetriggerthresholdof≥ 14hitPMTsarekept,resultinginthe neutrino effective area shown in Fig. 3. The neutrino effective area represents theequivalentdetectorareaforahypotheticalinstrumentwith100%efficiency for detecting the passage of neutrinos. The effective area is much smaller than the dimensions of the detector due to the small cross section of neutrino interaction and, for larger zenith angles and high energies, neutrino absorption on transit through the Earth. It is a useful parameter for determining event rates and making comparisons between experiments. The event rate for a neutrino model predicting a flux dΦ is given by dEνdΩν dN (cid:90) (cid:90) dΦ µ = dE dΩ Aeff(E ,Ω ) . (1) dt ν ν ν ν ν dE dΩ ν ν Using the atmospheric neutrino flux of Barr et al. [22], we find approximately 134 000 atmospheric neutrino events per year at trigger level with a median energy of 670 GeV and a slight density dependence on zenith angle, and thus declination, of ±15%. This number of events will be reduced however by topo- logical and reconstruction cuts necessary to eliminate the large background of misreconstructed downgoing muons from cosmic ray air showers. We esti- mate 50% of atmospheric events will be lost to reject this background; thus we choose 67 000 events from the atmospheric neutrino sample for our data set, roughly representing a year of data. Cuts applied to a real data sample may have some small energy dependence, i.e. efficiency may drop for low energy events, but we do not consider this effect. We generate astrophysical signal neutrinos similarly; however, they are simulated at a fixed declination and weighted differently than atmospheric neutrinos. We consider for the signal a 5 103 102 10 2 1 m 10-1 90o < Zenith < 100o 10-2 115o < Zenith < 125o 140o < Zenith < 150o 10-3 165o < Zenith < 175o 10-4 1 2 3 4 5 6 7 8 Log (E/GeV) 10 Fig. 3. Neutrino effective area at trigger level for several zenith bands. power law neutrino spectrum with spectral indices ranging from 1.5 to 4.0. The number of hit modules distribution for atmospheric neutrinos and several signal spectral indices is shown in Fig. 4. Angularreconstructionerrorsforacubickilometerdetectoriniceareshownto be∼0.7◦ forneutrinosofenergies(cid:38) 1TeV [1]. Wesimulatethisreconstruction error by adding to the muon direction a space angle error randomly sampled from a polar Gaussian distribution with standard deviation 0.7◦. Considering both this reconstruction error and the neutrino-muon vertex angle, the result- ing median angular resolution for an E−2 signal is 0.86◦ and is larger for softer spectra. A comparison to a detector with angular resolution 0.2◦, which may be achieved with a detector in sea water, is shown in Sec. 4. Muon energy reconstruction errors in neutrino telescopes are typically of the order of 0.3 in log E above a few TeV [23,24,25], limited by the stochastic 10 ν nature of muon energy losses. To simulate this effect, we assign to each event a reconstructed error sampled from a Gaussian distribution of width 0.3 in the logarithm of energy. The reconstructed energy for atmospheric neutrino events and several power law spectra is shown in Fig. 4. Energy resolution degrades below a few TeV; thus the energy resolution we assign to such events is unrealistically accurate. However, the power to detect astrophysical sources resides in the ability to discriminate high energy neutrinos from lower energy atmospheric neutrino background, so we conclude overestimating the energy resolution of such lower energy events does not affect the result significantly. 6 10-1 1 E-2 Neutrino E-2.5 Neutrino E-3 Neutrino 10-1 10-2 Barr et al. Atmos. Neutrino 10-2 d) E)] Nmo10-3 (og1010-3 dP/d(10-4 dP/d[l 10-4 E-2 Neutrino 10-5 E-2.5 Neutrino E-3 Neutrino 10-5 10-6 Barr et al. Atmos. Neutrino 10-7 0 20 40 60 80 100120140160180200 1 2 3 4 5 6 7 Number of Hit Modules (Nmod) Log (Estimated Energy/GeV) 10 Fig. 4. Number of hit modules distribution (left) and reconstructed muon energy distribution (right) for several neutrino spectra. 3 The unbinned likelihood ratio method Inthecontextofamuonneutrinopointsourcesearch,thedatafromaneutrino telescope consists of a set of muon events spread over the sky, each with reconstructed direction (declination and right ascension), energy, and time. While event time is useful in searches for short neutrino bursts or periodic neutrino emission, we focus on searches for continuous neutrino emission and disregard event time. The vast majority of events are muons produced by atmospheric neutrinos. At any celestial direction, the data can be modeled by two hypotheses: • H : The data consists solely of background atmospheric neutrino events. 0 • H : The data consists of atmospheric neutrino events as well as astrophys- S ical neutrino events produced by a source with some strength and energy spectrum. The likelihood of obtaining the data given each hypothesis is calculable, and the ratio of likelihoods, or equivalently the log of the likelihood ratio, serves as a powerful test. We define our test statistic (cid:34) (cid:35) P(Data|H ) 0 λ = −2·log . (2) P(Data|H ) S Largervaluesofλindicatethedataislesscompatiblewiththebackgroundhy- pothesis H . The probability density functions P(Data|H ) and P(Data|H ) 0 0 S are calculated using knowledge of the spatial and energy distribution of back- ground and astrophysical neutrino events. 7 Suppose we wish to test the existence of a source at a known direction(cid:126)x using s a set of data events, each with reconstructed direction (cid:126)x and reconstructed i energy E . Events reconstructed outside a declination band centered at (cid:126)x i s with a width several times the detector resolution are unlikely to have been produced by a source at (cid:126)x and can be disregarded, leaving N events in the s band. Each event inside the band is assigned a source probability density corresponding to the probability of the event belonging to a source at (cid:126)x . The s source is assumed to emit neutrinos according to an E−γ power law energy spectrum. The source probability density is the product of a spatial density function describing the potential of an event reconstructed with direction (cid:126)x i to have true direction (cid:126)x and the probability of observing reconstructed muon s energy E given source spectral index γ: i (cid:90) S ((cid:126)x ,(cid:126)x ,E ,γ) = N((cid:126)x |(cid:126)x )· P(E |E )P(E |γ)dE . (3) i i s i i s i ν ν ν Eν The spatial probability density component N((cid:126)x |(cid:126)x ) can be obtained directly i s fromeventreconstructionwhenmaximumlikelihoodreconstructiontechniques are used [26]. The structure of the reconstruction likelihood near the most likely direction (cid:126)x provides an event by event estimate of reconstruction un- i certainty and yields a Gaussian spatial probability density profile: Ni((cid:126)xi|(cid:126)xs) = 2π1σ2e−|(cid:126)xi2−σ(cid:126)x2s|2, (4) where |(cid:126)x −(cid:126)x | is the space angle difference between source and reconstructed i s event directions, and σ is the reconstruction error estimate. The vertex angle between the neutrino and muon is neglected in N ((cid:126)x |(cid:126)x ). This angle is negli- i i s gible compared to reconstruction error for cubic kilometer neutrino telescopes in ice with resolution ∼0.7◦ and high energy threshold of ∼100 GeV. Addition of the vertex angle error is discussed in section 4 for cubic kilometer telescopes in sea water with resolution ∼0.2◦. The integral over E is precomputed from ν detectorMonteCarlo.TheresultingtablesofP(E |γ)for1.0< γ <4.0,shown i for several indices in Fig. 4, are computed in steps of 0.01 in γ, interpolated linearly, and stored for reference.1 The resulting source probability density is Si((cid:126)xi,(cid:126)xs,Ei,γ) = 2π1σ2e−|(cid:126)xi2−σ(cid:126)x2s|2P(Ei|γ) (5) and has value unity when integrated over solid angle and E . The slight decli- i nation dependence of the atmospheric neutrino background can be neglected 1 For maximum accuracy, one may wish to tabulate P(E |γ) with respect to zenith i angle as well; however, we omit this step and achieve good results. 8 over the width of the band. The background probability density again depends on event energy and is then P(E |φ ) i atm B = . (6) i Ω band The probability density P(E |φ ) is precalculated as above assuming the at- i atm mosphericneutrinofluxofBarretal.[22].Weassumethebackgroundispurely atmospheric neutrinos. If the background contains another component, for ex- ample high energy muons from cosmic ray air showers, the energy component of the background density should be modified accordingly: P(E |φ +φ ). i atm Bkgd The source and background densities are combined, and the likelihood is eval- uated over all events in the band: (cid:32) (cid:33) n n (cid:89) s s L((cid:126)x ,n ,γ) = S +(1− )B (7) s s i i N N N where n describes the number of signal events present in the band. s The fraction of signal events n as well as the source spectral index γ are s not known and must be determined by maximizing the likelihood L. This is done by minimizing the quantity −log(L) using the MIGRAD minimizer available in MINUIT [27] with respect to the unknown quantities n and γ s and obtaining the best value of each parameter, nˆ and γˆ. The minimization s procedure also incorporates a penalty factor2 for γ > 2.7 to better discrim- inate astrophysical sources with hard spectral indices from the atmospheric background with γ ∼3.7. The original hypotheses can be written in terms of L: P(Data|H ) = L((cid:126)x ,0) and P(Data|H ) = L((cid:126)x ,nˆ ,γˆ). The test statistic 0 s S s s is (cid:34) (cid:35) L((cid:126)x ,0) s λ = −2·sign(nˆ )·log . (8) s L((cid:126)x ,nˆ ,γˆ) s s Ignoringthefactorsign(nˆ ),theteststatisticλisnevernegativesinceL((cid:126)x ,0) s s is contained in the range of L((cid:126)x ,n ,γ), of which L((cid:126)x ,nˆ ,γˆ) is the maximum. s s s s A downward fluctuation of the background may occur at (cid:126)x which would be s fit best as a source with negative number of events and a negative value of nˆ . s Such a downward fluctuation also would have a large value of λ, so sign(nˆ ) is s used to separate negative and positive excesses. A full sky search is a simple extension of this single point search method and can be accomplished by 2 The extension of γ is limited above 2.7 during the minimization by a Gaus- sian likelihood penalty with σ=0.2 in spectral index. This improves the discovery potential for hard astrophysical source spectra while remaining comparable to an unbinned method without energy for source spectra as soft as ∼3.2 (see Fig. 8). 9 performing the search on a grid of locations covering the sky. Finally, while it ispreferabletouseanenergyestimationtomaximizethepowertodiscriminate astrophysical neutrinos from the background, it is possible to do the search without this information. The energy dependent terms are removed (i.e., set to one) from the signal and background probability densities, Si((cid:126)xi,(cid:126)xs) = 2π1σ2e−|(cid:126)xi2−σ(cid:126)x2s|2 (9) 1 B = (10) i Ω band resultinginsimplerexpressionsforL.Thequantity−log(L)isminimizedwith respect to n , and P(Data|H ) = L((cid:126)x ,0) and P(Data|H ) = L((cid:126)x ,nˆ ). s 0 s S s s 4 Results We apply the likelihood ratio method to data consisting of the 67 000 back- ground events described in section 2 and an added source at declination 48◦. A grid of simulated source strengths and spectral indices is used, with source strength spanning 0 - 100 signal events added to the sample and spectral index γ spanning 1.0 - 3.9 in increments of 0.1, resulting in a total of 3000 simulated combinations of source strength and spectral index. For each combination, 10 000 experiments are done and the value of λ is recorded for each. 107 trials are performed with background alone to evaluate the significance of observed values of λ. The method is compared against a binned search with a circular bin centered at the source location. The bin radius is optimized by minimizing the number of E−2 signal events necessary to achieve 5σ significance in 90% of experiments. The optimal radius is found to be 1.35◦, with a signal efficiency of 80% and background expectation of ∼16 events/bin. Also, the method is compared against itself using only spatial information and neglecting event energy. An identical grid of simulated source strengths and spectral indices is used for both the binned method and likelihood method without energy. Fig. 5 shows a significance sky map with an added source at declination δ=48◦ and right ascension α=12h producing 15 events according to an E−2 energy spectrum. Fig. 6 illustrates the procedure used to determine significance and discovery potential. The integral distribution of λ for background alone is produced at declination δ=48◦ and the values of λ corresponding to 3σ (2.7× 10−3) and 5σ (5.7×10−7) integral probability are calculated. Fig. 6 also shows distributions of lambda with 8, 16, and 24 signal events added to background. Discovery potential at 5σ is the fraction of experiments with λ exceeding the 5σ threshold. Discovery potential is then computed in this fashion for each 10