Draftversion January10,2013 PreprinttypesetusingLATEXstyleemulateapjv.5/2/11 TIME-AVERAGE BASED METHODS FOR MULTI-ANGULAR SCALE ANALYSIS OF COSMIC-RAY DATA R. Iuppa1,2 & G. Di Sciascio2 Draft version January 10, 2013 ABSTRACT In the last decade, a number of experiments dealt with the problem of measuring the arrival di- rection distribution of cosmic rays, looking for information on the propagation mechanisms and the identification of their sources. Any deviation from the isotropy may be regarded to as a signature of 3 unforeseen or unknown phenomena, mostly if well localized in the sky and occurring at low rigidity. 1 ◦ It induced experimenters to search for excesses down to angular scale as narrow as 10 , disclosing 0 2 the issue of properly filtering contributions from wider structures. A solution commonly envisagedin these years based on time-averagemethods to determine the reference value of cosmic ray flux. Such n techniques are nearly insensitive to signals wider than the time-window in use, thus allowing to focus a the analysis on medium- and small-scale signals. Nonetheless, often the signal cannot be excluded in J thecalculationofthereferencevalue,whatinducesystematicerrors. Theuseoftime-averagemethods 9 recentlybroughttoimportantdiscoveriesaboutthemedium-scalecosmicrayanisotropy,presentboth in the northern and southern hemisphere. It is known that the excess (or the deficit) is observed as ] M less intense than in reality and that fake deficit zones are rendered around true excesses, because of the absolute lack of knowledge a-priori of which signal is true and which is not. This work is an I attempt to critically reviewthe use of time average-basedmethods forobservingextended features in . h the cosmic-ray arrival distribution pattern. p - o INTRODUCTION Amenomori et al.2004,2007;Abbasi et al.2011). Inthis r last case, as the amplitude and the phase of the signal t Large field of view (fov) experiments operated for s areanalyticallypredictable,theobservationiscommonly a cosmic-ray (cr) physics collect huge amount of high- considered as the starting point of any anisotropy anal- [ qualitydata,makingpossibletostudythecrarrivaldis- tribution with remarkable detail. Either satellite-borne ysis, as it demonstrates the reliability of the detector 1 and the analysis methods to be fine-tuned. Concerning or ground-based detectors are considered, many collab- v the “galactic” Compton-Getting effect, the importance orations coped with the measurement of the cr inten- 3 of this measurement lays in determining whether the cr sity all over the portion of the sky they observed. They 3 plasmaisco-movingornotwiththe“stillsystem”under all looked for deviations from the isotropic distribution, 8 study (Amenomori et al. 2006). as any signature of anisotropy provides essential infor- 1 Movingtonarrowerscales,itisknownthatacr“pure” mation on crs and the medium that they propagate 1. through. anisotropy,i.e. notduetoexpectedCompton-Gettingef- ◦ 0 Apart from the search for gamma-ray emission from fects, exists down to angular scales as wide as ∼60 . It 3 point-like(orquasi point-like)sources,eitherintheMeV wasobservedbyground-baseddetectorseversince1930s 1 energy range on-board satellites and in the TeV region and most recent experiments represented it in 2D sky- : maps (see Iuppa 2012 and references therein). Such a v with ground-based telescopes, directional data are ana- “large-scale” anisotropy (lsa) is of fundamental impor- i lyzed to map the cr gradient all over the sky at every X tance, commonly interpreted as a signature of the prop- angular scale. agation of crs in the local medium. r a Signal as deviation from the isotropy Beingchargedparticles,crshavetrajectoriesdeflected Any motion of the laboratory system with respect to by magnetic fields, so that their rigidity sets up a “mag- the cr plasma turns to a dipolar signature with a netic horizon”, i.e. a distance below which the ob- maximum in the direction of the motion. This is served arrival distribution contains information about true for any “still system” we want to consider: it the interaction of crs with the medium that they prop- might be well the Solar System (i.e. the motion agate through. The diffusion approximation effectively of the Earth around the Sun is factorized) or the explains the observations beyond this horizon. Thus, Galaxy itself (the motion of the Solar System around GeV-TeV crs are an effective tool to probe magnetic the Galaxy center is factorized). Such a process of fields within the Solar System (up to the Heliotail) dipole generationis commonly referredto as“Compton- (Desiati and Lazarian 2012). The multi-TeV region is Getting”effect(Compton and Getting1935)andwasob- important to study the cr propagation in the Local servedbya number ofexperiments (Aglietta et al.1996; Inter-StellarMedium(LISM),whereashigherenergycrs may reveal important features of the galactic magnetic [email protected] fields. If electrons (e±) are considered, synchrotron en- [email protected] ergy losses should be accounted for in defining the mag- 1Dipartimento di Fisica dell’Universit`a “Tor Vergata” di netic horizon (remarkably closer than for protons of the Roma,viadellaRicercaScientifica1,00133Roma,Italy. 2Istituto Nazionale di Fisica Nucleare, Sezione di Roma Tor same rigidity). Apart from that, the line is the same Vergata, viadellaRicercaScientifica1,00133Roma,Italy. and gives the importance of any attempt to measure 2 Iuppa et al. the anisotropy even in the e± channel. The cr elec- 2004). As well pointed out by (Fleysher et al. 2004), tron anisotropy was recently searched for by the Fermi these symmetries are assumed to be valid in some con- experiment (Ackermann et al. 2010), though with null ditions and such an assumption is part of the null- result. hypothesis against which signals are tested. Since 2009 “medium-scale” anisotropy (msa) struc- Among all the techniques that experimenters devel- tureswereobservedinthecrdistributiondowntoangu- oped to estimate the exposure, this paper is focused on larscaleaswideas10◦,theiroriginstillunexplained(cr thosebasedonthetime-averageofcollecteddata. What- source region, magnetic structures focusing crs, etc...) ever its particular implementation is, any time-average (Amenomori et al. 2007; Abdo et al. 2008a; Iuppa et al. method (tam) relies on the assumption that the signal 2012a). They were observed both in the Northern and to be interpreted as background can be filtered out by the Southernhemisphere atTeVenergyandthey donot averaging the event rate in a certain time window. The appear correlated to each other. It is quite natural to time-averagecan be performed directly, i.e. by integrat- expect that a well defined cr pattern at any angular ing the event rate within the time-window and then di- scale will be found as soon as statistics will be enough. vidingbythe windowwidth(direct integration method), Thenarrowerthe structurestheclosertheirorigin,what or by using Monte Carlo techniques. In the latest case, explains the growing interest towards this phenomenon eacheventisassociatedwithanumberN of“fake”events (Desiati and Lazarian2012; Drury and Aharonian 2008; havingall the same experimentalcharacteristicsbut dif- Salvati and Sacco 2008; Giacinti and Sigl 2012). ferent arrival time, sampled according to the measured Besidesthesecrsignals,diffusegamma-rayemissionis trigger rate. After “time swapping” (or “shuffling” or oftenmeasuredbysatelliteexperimentsintheGeVrange “scrambling”), the over-sampling factor N is accounted (Ackermann et al.2012)andextensiveair-shower(EAS) forandthefinalresultmakesthisapproachequivalentto arrays observed the diffuse emission from the Galactic the DirectIntegrationmethod. Actually, the onlydiffer- Plane at TeV energy (Abdo et al. 2008b). Structures as enceisthattheintegrationisperformedviaMonteCarlo ◦ ◦ wide as 10 − 20 have to be properly extracted from instead of directly. the background (the overwhelming cr contribution or Theuseof tamsisfavoredbythepropertywellknown the average photon content) and the analysis methods insignalprocessingforwhichsmoothingoutasignalwith usedto do thatoftenrestonthe same ideasexploitedto a top-hat kernel as wide as T strongly suppress signals measure the cr anisotropy. narrowerthanT. Asaconsequence,thesmoothedsignal will contain only signals wider than T, so that subtract- Detection techniques ing it from the actual signal is the same as saving only Either the cr anisotropy or the diffuse emission from contributions from frequencies higher than 1/T. As it the Galactic plane is considered, the experimental issue will be discussed in a more formal way in the next sec- of properly detecting and estimating the intensity of a tions, tams demonstrated to effectively work for point- signalasbrightas10−4−10−3withrespecttotheaverage like and quasi point-like gamma-ray sources, because in isotropic flux of crs has to be dealt with. these cases a suitable time interval around the source Ittranslatesinestimatingtheexposureofthedetector can be excluded from the integration. The time average withaccuracywellbelowthatthreshold,toavoidthatde- alongside the excluded region is quite a robust estima- tector effects mimic signals due to physics. To keep the tion ofthe exposure,i.e. - througha simple scalingwith exposure map under control down to ∼ 10−4 is a chal- the averagetrigger rate - of the cr background. lenge even for the most stable experiment, as unavoid- The attention on tams recently grew up because of able changes in the operating conditions occur and not their application to detect medium scale excesses on always on line corrections can be readily applied. Both top of the large scale cr anisotropy. The Milagro col- for satellite experiments and ground-baseddetectors the laboration firstly tried to “adapt” the Direct Integra- envisagedsolutionistoestimateoff-linetheexposure,re- tion method for studies on small to intermediate scales lyingonstatisticalmethods appliedtothelargedata-set (10◦-30◦) (Abdo et al. 2008a). The attempt came from available. the assertion that averaging for a time interval T cor- Good results can be achieved by combining Monte responds to make the analysis insensitive to structures Carlo simulations and the record of the operating con- wider than 15◦T/1 hr in Right Ascension. Afterwards, ditions (see for instance (Ackermann et al. 2012)). Oth- other experiments applied tams for anisotropy studies, erwise, if simulations do not reach the needed accuracy, either to estimate the over all exposure or to focus the like for EAS experiments3, or just to have a simulation- analysisonacertainangularscale(Guillian et al. 2007; independent result, the exposure is estimated from data Ackermann et al. 2010; Abbasi et al. 2011; Iuppa et al. themselves,by exploitingsome (assumed)symmetries in 2012a). Insomecases,thepropertyoffilteringoutlarger the data acquisition. structuresbecamethemainreasonwhy tamswereused, A number of data-driven methods to estimate the ex- although no detailed discussion was ever made on the posureexist,althoughallofthemarebasedongeometri- potentialbiasesofthese techniquesin filteringthe large- calpropertiesofthedetectoracceptance(seeforinstance scale structures. the “equi-zenith” methods (Amenomori et al. 2005)) This paper collects a series of simple calculations and and/or on the uniformity of the trigger rate within a observationsonthefilteringpropertiesof tams. Asthey certain period (Alexandreas et al. 1993; Fleysher et al. are applied in a variety of experiments having different operatingmodes, sky-coverageandtrigger rate stability, 3Theimportanteffectoftemperatureandpressurevariationson a generaltreatmentofthe matter is impossible. When a theatmosphericdepth and,consequently, onthetriggerefficiency specificexperimentallayouthadtobeaccountedfor,the ofEASarrayscannotbetakenintoaccountdownto10−4-10−3 authors made use of a virtual EAS array similar to the Time-average methods 3 ARGO-YBJ experiment (Bartoli et al. 2011a), whereas where N is the actual (unknown) background cr num- b todiscussalikelycaseofunderlyinglargescalestructure ber,N istheestimatedone,andN indicatesthenum- b ev the model of the large scale anisotropy of crs as given ber of measured events. The average is computed in the in (Amenomori et al. 2007) has been used. time interval T andusing the kernelfunction w, so that: e Statistical effects were not considered, i.e. no Pois- seoanchiapnixfleulcwtueraetiaocncsouanroteudndfort.heInafvaecrta,gtehiesvceonnttcriobnutteinotnoisf dNev(Ω,t) = tt−+TT//22dτ dNevd(τΩ,τ)w(τ) (2) tooutlinesomemajorpotentialsystematiceffects,intrin- (cid:28) dt (cid:29)w,T R tt−+TT//22dτ w(τ) sic to the application of tams, regardless if the number of events is sufficient to make them visible or not. If the source contribution is nRot excluded, the function The paper is organized as follows. In the section 1 w(τ) in (2) is the trigger rate and accounts for over-all an introduction on tams as exposure estimation meth- variations in the acquisition regime: ods is given. The section 2 is a brief interlude which dN (Ω,t) demonstratesaconsequenceofdata-normalizationalong dτ w(τ)=dτ dΩ ev (3) therightascensiontobeconsideredforallfurtherdiscus- (cid:20)ZFOV dt (cid:21)t=τ sions,thoughmostlyaffectingthe ℓ=1,2,3components The integration is carried out numerically, with the Di- of the signal. In the section 3 the effect of tams on the rect Integration or the Time Swapping method (see the signaltobedetectedareintroduced,mostlyforwhatcon- Introduction). cerns the reduction of the intensity and the appearance Thenextsectionsofthispaperwillfocusontheroleof of border effects. The section 4 finally provides quanti- dN /dτ in the estimate (2), as this quantity is the sum tative information on the residual contribution from fil- ev of different contributions and the problem of a proper tered components and the signal distortion due to the separation of the signal in the angular domain via tams method. Some conclusive remarks are given in the last has to be approachedby considering the time properties section. of dN /dτ. However, before that, two other aspects of ev the equations (1)-(2) should be made explicit. 1. EXPOSURECALCULATIONWITHTAMS From the experimental viewpoint, the observation of • Time interval. The quality of the approximation excess (or deficit) effects at a level of 10−4 is a difficult is related to the difference between N (Ω,t) and b task, because of the intrinsic uncertainty that cr appa- N (Ω,t) (1), i.e. to how representative the time- b ratus have to cope with in estimating the exposure. For averageis of each instant cr-flux. If the time win- EAS arrays the atmosphere is part of the detector it- dowT ischosentoowide,thegeometricaldistribu- e self and data must be handled with care to avoid that tion of the cr arrival directions may significantly a atmospheric change mimics a signal somewhere in the change, due to atmospheric effects. Some changes sky. For detectors on board satellite, no atmosphere ef- inthedetectoroperatingregimemayhavethesame fects are there but trigger rate variations persist related effect of making the dN /dτ distribution not uni- ev to changing conditions along the orbit. form. In general, assuming there is an isotropic charged cr flux overwhelming all the other signals, the exposure is • Source exclusion. The source contribution should estimated by assuming it proportional to the integrated be excluded from the time-average. Mathemati- cr flux. In this way,the exposure estimation problemis cally, the weight (3) has to be replaced by w (τ): se posed as a cr-counting problem. Hereafter, the number of events collected (or com- w(τ)−→w (τ)= 0 if Ω∈Dsrc(τ) (4) puted) in the solid angle dΩ centered around Ω =(θ,φ) se w(τ) otherwise (cid:26) in the local frame, in the time interval [t,t+dt) will be written as: where Dsrc(τ) indicates a confidence solid angle dN(Ω,t) d2N(Ω,t) around the source at the time τ. = dt dΩdt As far as the time window is concerned, the acquisi- (cid:18) (cid:19) tion of eas arrays is not stable for periods longer than to lighten the notation. 2−3 hrs,as climatic changesaffect either the arrivaldi- rection distribution of cr and the detector response to 1.1. Point-like and quasi-point-like sources the incoming radiation. There are far minor problems Forpoint-likeorquasi-point-likesources,tamsareusu- for underground experiments or neutrino observatories, whereevenlongertimes areusedin the literature(up to ally applied to estimate the exposure (i.e. the expected backgroundcr rate) from a certain direction of the sky. 24 hrs (Abbasi et al. 2011)). Nonetheless, time intervals as short as 4 hrs or less are used also in some of these They are an evolution of the elder “on-off” method and rely on the assumption that the cr flux froma givendi- cases to extract small-scale signals. Satellite-borne de- tectorsusuallyaresostabletoallowtotheexperimenter rection Ω in the local reference frame is practically con- to shuffle events within the whole data-set available (up stant during short time-periods. In other words, the av- to few years), so that the analysis dows not suffer the erage count from Ω =(θ,φ) during the the interval T is quite a good approximation of the cr number N : pitfalls described below (Ackermann et al. 2010). cr About the source exclusion, the solid angle to be ex- dN (Ω,t) dN (Ω,t) dN (Ω,t) cluded around the source is related to the detector an- b b ev ≃ = (1) gularresolution. A safe choice might be 2 or3 times the dt dt dt (cid:28) (cid:29)w,T e 4 Iuppa et al. average angular resolution plus the source intrinsic ex- extension T in local hour angle and time-window T it S tension. If a 2◦-wide source is observed with an angular holds ρ=1−T /T. S ◦ ◦ ◦ resolutionof1 ,asafeexclusionregionof6 −8 around it can be set. If the region is populated of other known 2. TAMSANDLSA sources, the definition of the exclusion region has to be Before coping with the filtering properties of tams, obviously adapted. we discuss here the effect of tams on the measurement For all experiments surveying the sky, the fov does of the lsa of cosmic rays. Actually, no modern ex- not coincide with the portion of the celestial sphere to periment but IceCube used tams to estimate the expo- be investigated. They exploit the rotationof the labora- sure (Abbasi et al. 2011) for all-scale analysis, because tory frame with respect to the sidereal frame to get the it would mean to average along 24 hrs and to face all project coverage. In this sense, all time-spans may be the issues of detector stability addressed in the previous translated into angular intervals measured in the side- section. Nonetheless, the result reported here is valid real frame. If the laboratory rotates around the Earth alsoforallthe othermeasurementsofthecr anisotropy, axis (ground-basedexperiments), time intervals arer.a. e.g. “equi-zenith” (Amenomori et al. 2005) or “forward- intervals. Depending on the rotation of the laboratory backward”(Abdo et al. 2009)orelse. Infact,acommon ◦ frameinthesiderealframe,1hrmaycorrespondto∼15 device to bypass the ignorance of the absolute detection in r.a. for a ground-baseddetector or ∼240◦ in the or- efficiency as a function of the arrival zenith angle (i.e. bit panle of a low Earth circular orbit satellite. For the of the declination), is to set the average flux of cosmic IceCube detector,atthe SouthPole,the skyportionob- rays detected in a certain zenith (declination) belt to a served is always the same and time-flow simply brings certain value, the same for all different belts. In other a rotation with respect to the celestial coordinates. For words, deviations from the isotropy are not measured ◦ ◦ ground-baseddetectors2−3hrscorrespondto30 −45 withrespectto the averageoverthe whole skyobserved, andenoughstatisticsislefttoallowthe sourceexclusion astodothattheefficiencyofthedetectorasafunctionof (∼50−80% of the events inside the time window T can the zenith must be properly accounted for. Conversely, be used). In the literature, typical values are found to thereferenceaverageiscomputedalongeachzenithbelt. ◦ beT =2hrsforthetimeintervaland∆=6 fortheex- We show here that this solution introduces a degener- clusion region width (Fleysher et al. 2004; Bartoli et al. acy in the measurement of the anisotropy, i.e. m = 0 2011b). components of the signal are suppressed. Ifthesignalislookedatasadistributionf(θ,φ)onthe 1.2. Wider structures sphere,theactofnormalizingtheaveragecontentofeach Ifwiderstructuresareconsidered,thetwoconditionsof declination belt to zero can be written as the operator the previoussectioncannotbe fulfilledat the same time. S: In fact, the off-source integration interval becomes nar- 1 2π rowerthan the source extension, thus making the on/off f(θ,φ)−→S f′(θ,φ)=f(θ,φ)− dφf(θ,φ) source event ratio too high and introducing large fluctu- 2π Z0 ations in the exposure estimation. ′ where the f distribution is the measured one, which This is true for a number of structures having physi- differs from the “true” f for the average hfi = calextensionwiderthanfewdegrees. Forinstance,when θ 2π experimentslikeMilagroorARGO-YBJmeasurethedif- 1/2 dφf(θ,φ). We can consider the spherical har- 0 fuse emission from the Galactic Plane, the source exclu- monics expansion of the f distribution: sionregionisusuallya±5◦ galacticlatitude beltaround R ∞ ℓ the plane. Studies of systematics are performed by ex- ◦ f(θ,φ)= aℓ Ym(θ,φ) tendingtheregionupto±10 ,obtaininganon-negligible m ℓ contribution to the uncertainty (∼ 10% (Abdo et al. ℓ=0m=−ℓ X X ◦ 2008b)). The±5 choicegiveslessfluctuationsbutprob- andconsideringinacloserdetailtheeffectoftheaverage ably still includes some signal events in the background ◦ on the signal. In fact: estimation. On the contrary, the ±10 is a safer choice forwhatconcernsthesourceexclusion,atexpenseofthe 2π 0 if m6=0 statistics4. Quotingthis effect asa sourceofsystematics dφYm(θ,φ)= ℓ Ym(θ,φ) if m=0 is still acceptable because the experiments do not have Z0 (cid:26) ℓ ◦ ◦ thesensitivitytoextendthemeasurementupto10 −15 ′ Using the last result, f can be rewritten as: fromtheGalacticPlane. Perhapsnext-generationexper- iments will have it and it will not be possible to exclude ∞ the whole region of interest when applying tams. f′(θ,φ)=f(θ,φ)− aℓY0(θ,φ) (5) 0 ℓ AsimilarpointholdsfortheMSAregions,oftenwider ℓ=0 ◦ X than20 ,forwhichthesourceexclusionisnotapplicable. where the degeneracy is made explicit. In fact, allterms In these cases, the signalintensity is reduced by a fac- withm=0aresuppressedbytheexperimentaltechnique tor ρ depending on the signal and the background mor- applied, what is more important as the multipole order phology, as well as on the time-window chosen to apply the tam. For uniform background, uniform source with ℓ gets lower. Iftheskyisonlypartiallyobserved,furthereffectsarise 4 It should be noticed that ±10◦ in Galactic Latitude corre- due to the non-uniform exposure. In fact, if the number sponds to a varying r.a. interval, as the Galactic plane is not of events strongly depends on the declination or other orientedalongthecelestialequator. preferred directions, significant deviations from isotropy Time-average methods 5 (a) (a) (b) (b) (c) (c) Figure 1. Effect of the r.a. normalization on a dipole signal. Figure 2. Effect of the r.a. normalization on a quadrupole sig- Figure (a) represents the input dipole signal, as intense as 0.1. nal. Figure(a)representstheinputquadrupolesignal. Figure(b) The dipole vector points at (θ = 63◦,φ = 243◦). Figure (b) rep- representsthequadrupolereconstructedwithmethodsnormalizing resentsthedipolereconstructedwithmethodsnormalizingtheav- theaveragecontentineachdec. belttozero. Figure(c)represents erage content ineach dec. belt to zero. The dipole vector points the difference between the input map and the reconstructed one, at(θ=90◦,φ=243◦)andtheintensityis0.089. Figure(c)repre- whichturnsouttobeproportionaltotheY02(θ,φ)function. sents the difference between the inputmapandthe reconstructed one,whichturnsouttobeadipoleasintenseas0.045,pointingat compute for each point the average: θ=0◦. Noticethat0.0892+0.0452=0.12. n+N/2 1 might be observed only in certain regions of the field of ξ = x (N ≤N) (6) n k view. N k=n−N/2 A representation of the effect just described is given X for ℓ=1,2 in figures 1-2. k=k±N if k <1 or k>N)thenthedifferencex−ξwill maintain intact all structures narrowerthan N, whereas 3. FILTERINGPROPERTIESOFTAMS all features much wider than N will be suppressed. In As time is a synonym for r.a., tams average signals Fig. 3 we show the results of a toy numericalestimation along the r.a. direction, i.e. they enjoy the property of of the time-averageeffect on the signal intensity estima- filtering out large scale contributions to the signal. tion for different angular scales. The red curve clearly If we have a data series x (k = 1,2,...,N) and we shows that the average along ∆T preserves signals on k 6 Iuppa et al. narrower angular scales and strongly reduces wider con- 1. theydonotinheritsystematicsfromeffectspresent tributions. The Fig. 3 triggers some other considera- below the angular scale they are set to filter out. tions. Firstly the excess (or the deficit) is observed as In this sense, no Compton-Getting interference is less intense than it really is. This bias can be avoided expected, neither influence or artifacts induced by by excluding the source region, what is impossible for large scale atmospheric effects, which instead were structureswiderthanhalfthetime-windowextension. A demonstrated to be relevant for the lsa analysis; second important issue (related to the first one) is that In fact, systematics introduced in misinterpreting fakedeficit zones arerenderedaroundtrue excessesand, the detector performance usually affects the whole vice-versa, fake excesses are seen around true deficits. sky,hardlybeingresponsibleforlocalizedfeatures; Theimportance ofthis problemlaysinthe absolutelack ofknowledgea-priori ofwhichsignalistrueandwhichis 2. the amplitude of the systematics described above, i.e. the residuals from the lsa structures, can be not. The problem is present mostly for structures com- evaluatedwithMonteCarlosimulationsifindepen- ing from reducing wider features than for really-narrow dent measurementsor robustmodels areavailable. signals. If the actual signal is well separated in the har- monic space from the wider structure to be suppressed, 4.1. Residuals of underlying large-scale structures it may be possible to observe it also without any fil- ter, i.e. just by considering the r.a. distribution like Beforeconsideringhowtamsreconstructextendedsig- in (Abdo et al. 2008a). On the other hand, if some hy- nals like the MSA, i.e. which effects the analysis tech- pothesis can be made (from the literature or from in- nique has on the intensity and the shape of the signal dependent data about the detector exposure), the effect under observation, it is worth to investigate the residual of such underlying larger scale signals can be easily esti- effect of the underlying large-scale structures, which are mated and quoted as systematic uncertainty of the final strongly suppressed by the time-average. measurement. Signal gradient along the RA direction 4. TAMSANDMSA The result of the tam depends on the signal which it is The study of the msa in the arrivaldistribution of cr applied to, so that no prediction is possible if the signal in toto isnotconsidered. Nonetheless,itisconvenientto can be approached with different methods. focus on two characteristics of the signal separately, i.e. Likely the most orthodox way is to evaluate the expo- the extension and the gradient along the RA direction. sure with techniques sensitive to any angular scale and It is easy to assess that the gradient of the estimated then apply the spherical harmonics analysis to filter out backgroundisproportionaltothedifferenceofthesignal the signal. The use of the a coefficients prevents any ℓm gradient at the boundaries of the time average window: contamination from harmonic regions other than those selectedwiththe(ℓ,m)numbersandallowstodefinethe d dN (Ω,t) ev degree of anisotropy in a (mathematically) robust way. ≃ dt dt Another approach, still starting from an all-scale- (cid:28) (cid:29)w,T sensitive estimation of the exposure, is to estimate the k dN (Ω,τ) dN (Ω,τ) dipole and quadrupole components of the measured cr ev ev − (7) distribution,tosubtractthemfromtheexperimentaldis- T dτ (cid:12)τ=t+T/2 dτ (cid:12)τ=t−T/2! tributionandtofocusontheresidualsatscaleslessthan (cid:12) (cid:12) 90◦. For such narrower signals, the analysis is carried The constant k de(cid:12)(cid:12)pends on the kernel fun(cid:12)(cid:12)ction used5 In out in the real domain. the numerical implementation, if a “top-hat” kernel is These two methods enjoy the uncontested feature of used, it turns out to be k = 1. The equality does not properlyfilteringthesignalfromscaleslargerorequalto hold because of the denominator in the equation 2, in- 90◦, although some problems are there due to the par- troducingsecondordercorrectionsindNev(Ω,t)/dt. The tial coverage of the sphere by the experimental data. In equation 7 can be figured out by thinking of a discrete fact,eitherthea expansionandthedipole-quadrupole implementation of the tam, where a moving time win- ℓm determination are achieved with fit procedures over the dow passes from the i-th time-bin to the (i+1)-th. The whole sphere: that part of the sky which is empty of content of the bin centered at ti+1+T/2 is included in data has to be suitably masked and the lack of informa- the backgroundestimation in place of the contentof the tionunavoidablyreflectsupontheerrorassociatedtothe bin centered at ti−T/2. intensity and the structure position on the sphere. A simple representationof this effect is given in figure Moreover, both the a and the dipole-quadrupole 4. Thetoppanelrepresentalargescaleexcessasintense methods rely on an estimℓmation of the exposure all over as I =10−3 with respect to the isotropic cr flux, drawn the angular scales, what implies that they can filter out according to the equation: the lsa only if it is properly detected. If the all-scale dN (Ω,t) analysis revealed some systematics for the lsa, it would ev = be difficult to makesensefromana expansionofsuch dt ℓm a signal. I t−α +∆α/2 t−α −∆α/2 0 0 On the other hand, we already hinted that the appli- tanh −tanh 2 w w cation of time-average techniques to get the msa signal (cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21) would introduce systematics on the flux estimation, as The signal center is fixed at α = 12 hrs, the width 0 wellasthebiasoffilteringonly alongther.a. direction. at ∆α = 6 hrs and the signal gradient constant 1/w is Nevertheless,twoargumentsareinfavorofthesetech- niques: 5 Itmustbesymmetricintime,i.e. w(τ)=w(−τ). Time-average methods 7 · 10-3 s s e 0.1 c x e e v ti 0.08 a el r 0.06 0.04 D T 0.02 0 0 2 4 6 8 10 12 14 16 18 20 22 24 r.a. (hrs) (a) ·10-3 ·10-6 relative excess000...000567 relative excess 108 0.04 6 0.03 4 0.02 0.01 2 0 0 -0.01 14.2 14.4 14.6 14.8 15 15.2 15.4 15.6 15.8 19.6 19.8 20 20.2 20.4 20.6 r.a. (hrs) r.a. (hrs) (b) (c) Figure 3. Simulation of the effect of the time average on the estimate of the intensity of signals of different angular extensions. The central highestsignalpeaks around12hrs,another signalisvisibleinthespikearound15hrsandisrepresentedbythebluecurveinthe zoom (b). Onemoresignal,around20hrs,isrepresentedalsointhezoom (c). Black curve: sumofallthesignals. Red curve: the signal oftheblackcurveasitwouldappear aftertheapplicationofa3hrs-widetam. Blue curve: inputsignalcenteredatr.a.=15hrs. changed from ∞ to 1/1.5 hrs−1 (black to green curves). Signal gradient Amplitude of Width of (hrs−1) residuals(×10−4) residuals(hrs) The bottom panel reports how a tam (top-hat kernel, ∞ 5.0 1.5 time-window 3 hrs) would filter each signal of the panel 1/0.25 3.3 1.7 above. Residuals of the large scale signal remain, whose 1/0.5 2.2 2.2 intensity strongly depends on the signal gradient. For 1 1.0 3.0 a non-physical signal like the black one (w = 0), bor- 1/1.5 0.5 3.0 der effects as intense as half the input signal are visible. These effects are reduced below 10% the input intensity Table 1 if the signal gradient along r.a. is less than 1 hr−1. It Intensityandextensionofresidualsreportedinthebottom panel offigure4. can be noticed that residuals are both of positive and negative nature. The maximum intensity of the resid- The effect of the LSA on the time-average expo- uals and their extension (intended as the width of the sure estimation intervalswhereinthey arealwaysaboveorbelow 5%the Inthelastsectionitwasshownhowimportantisthegra- inputsignalintensity)arereportedinthe table1. Itcan dient of the large scale signal intended to be filtered out be noticed that the width is never above 3 hrs, i.e. the withtams. Resultsweregivenforatoy-calculationwith time-window used in the analysis. thepurposeofenhancingpotentialbiasesoftheanalysis. The aim of this section is to evaluate the effect of the 8 Iuppa et al. wasreplacedwiththeaveragecontentofallpixellessdis- 1.001 tantthanT/2. Thetam wasrepeatedforallvaluesofT 1.0008 from 1 to 24 hrs. Results for 1, 2, 3, 4, 6, 12 and 24 hrs 1.0006 will be reported only. e 1.0004 The Healpix “ring” pixelization scheme was used 1.0002 (Gorski et al. 2005). 1· 10 0· 10-32 4 6 8 10 12 14 16 18 20 22 24 Anticipatingthesection4.2,wenoticethattheangular 0.4 distance between pixels is the same all over the sphere, 0.2 (e-b)/b 0 bcountsitdheerart.ioan. adriestcalnocsee itnocrtehaesepsolwesh.enCtohnesepqiuxeenltulyn,daers -0.2 the interval for the average is set in the r.a. space, -0.4 the computation at low dec. values is carried out on 0 2 4 6 8 10 12 14 16 18 20 22 24 more pixels than at higher dec.. The methods becomes a [hrs] ineffective when the number of pixels at a certain dec. is such that the average along a certain ∆T is the pixel Figure 4. Toy-calculationtohighlighttherelevanceofthegradi- ent along the RA direction of the signal under study. Top panel: itself, making the background equal to the signal. The a6-hrswidesignalasintenseas10−3 withrespecttotheunderly- effect is small for dec. δ < 45◦ and T > 3 hrs = 60◦, ing flat background was simulated to riseup with different slope; but is important above δ =70◦, mostly for T ≤2 hrs. bottompanel: thesignalasreconstructedafterthetam-calculated Figure 5 reports the result of the calculus for different backgroundsubtraction(time-window: 3hrs). dec. bands. Theplotswereobtainedbyprojectingdata residualcontributionofalargescalesignalwhosenature ofthe eventmap(dN /dtabbreviatedwitheinthe fig- ev is closer to reality. Once again, the result was achieved ure) and the estimated background map (dN˜ /dt → b) with numerical calculations, and the algorithm applied b in the declination interval indicated, then by calculating is described in the following. for each bin the ratio (e−b)/b. In every plot, the black lsa parameterization.— To avoidto introduce anycircu- curve representsthe input lsa signalafter the r.a. nor- larbias,weusedtheparameterizationofthelsagivenby malization. Asnodetector-inducedeffectsnorstatistical theTibet-ASγ collaborationin(Amenomori et al.2010). fluctuations are considered, the gray curve perfectly fits Theauthorsmodeltheirobservationwithtwostructures, with it, representing the signal as it would be observed a Global and a Medium Anisotropy. The former signal with an all-scale sensitive tam (T = 24 hrs). The other is whatis commonlyreferredto as lsa,whereasthe sec- curvesrepresentwhatremainsofthelsastructurewhen ond one lays on smaller scales and is part of the signal the T-wide filter is applied. It can be appreciated how that the methods discussed in this paper are tuned for. the lsa is practically suppressed for T ≤4 hrs. The ta- It is worth noticing here that the best fit of the model ble2reportstheabsolutemaximumdeviationofthelsa- (Amenomori et al.2010),isgivenafterthenormalization inducedsignalfromthezero-reference-valueasafunction of data along each declination band (see the section 2). of the dec. band and the time-averagewindow, in units of10−4. Itisworthnoticingthatthemaximumeffectoc- Exposure simulation.— The detector exposure was sim- curs at the maxima of the lsa, so that for medium scale ulated by using the local cr distribution obtained for structures observedin other regionsthan those maxima, a flat ground-based detector, with a standard atmo- such an effect is far less than reported in the table 2. sphere absorption model (dN/dθ = I exp(−k/cosθ), 0 with k = 4.8). After time discretization, the local cr dec. Systematic distribution was computed for each time bin in the side- band error (×10−4) real day and transposed in equatorial coordinates. The 1hr 2hrs 3hrs 4hrs 6hrs exposuremapshowedthecharacteristicfeatureofamax- −20◦↔0◦ 0.07 0.29 0.7 1.1 2.4 imum at a dec. few degrees above the experiment lati- 0◦↔20◦ 0.06 0.26 0.6 1.0 2.1 tude (30◦ N), fading away at the field of view dec. lim- 20◦↔40◦ 0.05 0.20 0.4 0.8 1.6 its . As the calculus was performed with events arriving 40◦↔60◦ 0.04 0.14 0.3 0.5 1.1 withinθ =50◦,thereferencedec. bandforthisanalysis 60◦↔80◦ 0.016 0.06 0.14 0.25 0.5 ◦ ◦ is −20 ↔80 . Table 2 Systematicerrorinducedbythetime-averagefilteringmethod Event simulation.— As already said, in order to exclude usedintheanalysis. Thelsawasparameterizedaccordingtothe uncertainties due to statistics, the event map was not “Globalanisotropy”givenin(Amenomorietal.2010). filled by following a Poissonian distribution, but using the average value expected from the numerical integra- Figure 5 makes a point in showing how weak the ef- tionofthe localdistributionfunction. Thisis the reason fect of the lsa is when structures narrowerthan 45◦ are why no fluctuations are visible in the r.a. profile. The looked for. If the numbers of the table 2 are considered, actual background map wasobtainedwiththesamesolu- together with the typical intensity quoted for the msa tion: to avoid any fluctuations, any pixel was filled with emission (4· 10−4 −10−3 relative to cr isotropic flux, the product of the exposure times the total number of with T ≤ 3 hrs (Abdo et al. 2008a; Iuppa et al. 2012a; detected events. Abbasi et al. 2011)), systematic residuals from the lsa Estimatedbackground.— Theestimatedbackgroundmap after the tam are proved to be less than 20%. was obtained from the event map, i.e. from the sky pic- 4.2. Reconstruction of MSA structures turecontaining thelsasignal. Thecontentofeachpixel Time-average methods 9 (e-b)/b 0.00.002052 (e-b)/b 0.00.002052 0.0015 0.0015 0.001 0.001 0.0005 0.0005 0 0 -0.0005 r.a. projection d=-20(cid:176) - 0(cid:176) -0.0005 r.a. projection d=0(cid:176) - 20(cid:176) total MC input total MC input -0.001 1 hr -0.001 1 hr 2 hrs 2 hrs 3 hrs 3 hrs -0.0015 4 hrs -0.0015 4 hrs 6 hrs 6 hrs -0.002 12 hrs -0.002 12 hrs 24 hrs 24 hrs -0.0025 -0.0025 0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350 r.a.[(cid:176)] r.a.[(cid:176)] (a) (b) (e-b)/b 0.00.002052 (e-b)/b 0.00.002052 0.0015 0.0015 0.001 0.001 0.0005 0.0005 0 0 -0.0005 r.a. projection d=20(cid:176) - 40(cid:176) -0.0005 r.a. projection d=40(cid:176) - 60(cid:176) total MC input total MC input -0.001 1 hr -0.001 1 hr 2 hrs 2 hrs 3 hrs 3 hrs -0.0015 4 hrs -0.0015 4 hrs 6 hrs 6 hrs -0.002 12 hrs -0.002 12 hrs 24 hrs 24 hrs -0.0025 -0.0025 0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350 r.a.[(cid:176)] r.a.[(cid:176)] (c) (d) (e-b)/b 0.00.002052 0.0015 0.001 0.0005 0 -0.0005 r.a. projection d=60(cid:176) - 80(cid:176) total MC input -0.001 1 hr 2 hrs 3 hrs -0.0015 4 hrs 6 hrs -0.002 12 hrs 24 hrs -0.0025 0 50 100 150 200 250 300 350 r.a.[(cid:176)] (e) Figure 5. lsasimulationfordifferentdec. bands: (a)−20◦↔0◦;(b)0◦↔20◦;(c)20◦↔40◦;(d)40◦↔60◦;(e)60◦↔80◦. Different curvesrepresentdifferenttime-averagewindows(seethelegendfordetails). As the signal to be resolved is not excluded in the stant,whatisequivalenttoassumingthatthelargerscale background computation, distortions will appear in its structuresarewellseparatedintheharmonicspacefrom reconstruction. From previous section it should be clear the msa. In this sense, it should be noticed that the thattheeffectofthetamdependsontheactualshapeof calculation is the same of that reported in the study of the event distribution, whose composition is not known ther.a. gradient,withtheimportantdifferencethatthe a-priori, i.e. nobody knows what is signal and what is ratioof the signalwidth to the time window was greater background. Afterconsideringhowtheunderlyingwider than 1 there, but it is smaller here. structures affect the analysis, we describe here the re- The issue of the actual angular size of the signal and lation between the actual and the reconstructed signal. that of the reduction of the intensity are addressed. The results are obtained by considering a top-hat one- dimensionalsignaltowhichtamwithdifferenttimewin- Signal extension and declination dowsareapplied. Thebackgroundisassumedtobecon- If it is true that time is the same of RA, so that time- 10 Iuppa et al. average corresponds to RA-average, it does not mean be used has to be externally directed (i.e. after experi- that the filtering properties of tams are the same all mental trials or other measurements from literature). It overthe sphere. In fact, the angle ψ between two events is worth recalling that the relation 9 is obtained for a having coordinates (α ,δ) and (α ,δ) depends on δ: particular (and non-physical) case and it gives an upper 1 2 limit of the dependence of the signal reduction on the cosψ =sin2δ+cos2δ cos(α1−α2) (8) quantities ρ and α. The equation 8 clearly shows the dependence of ψ on δ CONCLUSIONS ◦ ◦ (as expected, ψ = 0 if δ = 90 ). As a consequence, This paper is intended to contribute to the study of the effect of any filter working in the RA space will be the analysis methods implemented to observe extended different according to the declination band: a top-hat ◦ ◦ signals in the cosmic ray arrival direction distribution. filter45 wideinRA,correspondstoatop-hatfilter31.4 Inparticular,itfocusesonthe time-averagebasedmeth- wide. A representation of the equation 8 for α −α = ◦ 1 2 ods applied when the source cannot be excluded, and 45 is given in figure 6. it points out how they undoubtedly enjoy properties of filtering in the r.a. direction, at expense of important spurious effects introduced in the signal reconstruction, ] 50 (cid:176)[ both for what concerns its intensity and its shape. Ex- y 45 perimentslikeMilagro,IceCubeandargo-ybjmadeuse 40 of these methods in the last decade (Abdo et al. 2008a; 35 Iuppa et al. 2012a; Abbasi et al. 2011). On the important point of potential residual effects 30 coming from larger signals, supposed to be filtered out, 25 it has been shown that they depend on the gradient of 20 such signals, rather than on their intensity, what makes 15 the bias from the known underlying lsa really small. Knowndistortionsofthereconstructedsignalwerean- 10 alyzed,givingnumericalinformationaboutthe intensity 5 reduction,the presence of bordereffects (deficits around 0 -80 -60 -40 -20 0 20 40 60 80 excessesandvice-versa)andasortof“reshaping”dueto d[(cid:176)] the filter acting only along the r.a. direction. Figure 6. Angular distance between two points on the sphere We conclude that the detection of medium and small having RA coordinates shifted45◦ fromeach other, as a function structures with tams may hardly mimic fake signals ofthedeclination. Theverticallinesenclosethedeclinationregion more intense than 210−4 with respect to the average forwhichmanycomputations ofthispaperweremade. isotropic cr flux. Nonetheless, if more detailed studies are attempted, like energy spectrum (i.e. signal inten- This is the reasonwhy tams tend to return structures sity) or morphology, care is needed in considering sys- moreandmorenarrowasdeclinationbandsfartherfrom ◦ tematic effects introduced by tams. In the near future, δ =0 are considered. experiments(Abbasi et al.2011;Abeysekara et al.2012; Cao et al.2010) willhavethe sensitivityto go below the Signal reduction We consider a signalwith intensity dN /dt and width w 10−4 level, what will make the choice of tams less and s (inr.a.),aboveaflatbackgroundwithintensitydN /dt. less effective. b It is analyzed with a tam with time window T. The In this sense, medium and small-scale cr anisotropy (measured)eventcontentis indicated with dN /dt. The are best searched for if a full-scale analysis is applied, e reconstructed signal a′ withstandardsphericalharmonicsorwavelettechniques (Iuppa et al. 2012b), allowing the experimenter to focus ′ dNe/dt − dNb′/dt on well-defined regions in the harmonic space. The pos- a = dNb′/dt sibility of this approach is related to the capability of controlling the detector exposure down to the level of can be compared with the actual one the signal to be observed, i.e. accounting for detector andenvironmenteffectsaspreciselyas10−4−10−3. Due dN /dt − dN /dt dN /dt a= e b = s to this reason, time-average methods, if carefully imple- dNb/dt dNb/dt mented, can still be considered a good compromise for ′ current experiments, which are controlled to this level by studying the ratio α = a/a. The ratio will depend but do not have yet the sensitivity to detect signals as on the ratio ρ = w/T and on a itself. It is easy to low as 10−5. demonstrate that for 0 < ρ ≤ 1 and ρa ≪ 1 (conditions fulfilled for the measured msa intensity) it holds: REFERENCES α≃(1−ρ)(1+ρa)≃1−ρ (9) R.Abbasietal.,Apj,740,16(2011). i.e. the reductionofthe signalintensitydepends linearly A.A.Abdoetal.,Physical Review Letters,101,221101(2008). on the ratio ρ. A.A.Abdoetal.,ApJ,688,1078(2008). A.A.Abdoetal.,ApJ,698,2121-2130(2009). The equation9 makesexplicit one of the majoruncer- A.U.Abeysekaraetal.,Astroparticle Physics,35(10), 641-650 tainty induced by tams: since the signal width is some- (2012). thingtobedetermined,thechoiceofthetimewindowto M.Ackermannetal.,Physical Review D,82,092003(2010).