A New Maximum-Likelihood Technique for Reconstructing Cosmic-Ray Anisotropy at All Angular Scales

3.5 MB·

by  Markus Ahlers

PreprinttypesetusingLATEXstyleemulateapjv.08/22/09
A NEW MAXIMUM-LIKELIHOOD TECHNIQUE FOR
RECONSTRUCTING COSMIC-RAY ANISOTROPY AT ALL ANGULAR SCALES
M.Ahlers1,†, S.Y.BenZvi2, P.Desiati1, J.C.D´ıaz–Ve´lez1,3, D.W.Fiorino1,4, and S.Westerhoff1
ABSTRACT
The arrival directions of TeV–PeV cosmic rays show weak but significant anisotropies with rela-
tive intensities at the level of one per mille. Due to the smallness of the anisotropies, quantitative
studies require careful disentanglement of detector effects from the observation. We discuss an itera-
6
tive maximum-likelihood reconstruction that simultaneously fits cosmic ray anisotropies and detector
1
acceptance. The method does not rely on detector simulations and provides an optimal anisotropy
0
reconstructionforground-basedcosmicrayobservatorieslocatedinthemiddlelatitudes. Itisparticu-
2
larlywellsuitedtotherecoveryofthedipoleanisotropy, whichisacrucialobservableforthestudyof
y cosmicraydiffusioninourGalaxy. Wealsoprovidegeneralanalysismethodsforrecoveringlarge-and
a
small-scale anisotropies that take into account systematic effects of the observation by ground-based
M
detectors.
Subject headings: cosmic rays — reference systems — methods: data analysis
8
1
1. INTRODUCTION Mertsch 2015; L´opez-Barquero et al. 2015). The possi-
] ble influence of the heliosphere via magnetic reconnec-
E During the past decade, a number of cosmic-ray, γ-
tions in the heliotail (Lazarian & Desiati 2010), non-
H ray, and neutrino observatories have found anisotropies
isotropic particle transport in the heliosheath (Desiati
in the arrival directions of Galactic cosmic rays at TeV
. & Lazarian 2013) or the heliospheric electric field struc-
h andPeVenergies(Guillianetal.2007;Amenomorietal.
ture(Drury2013)hasalsobeenconsidered. Moreexotic
p 2006; Abdo et al. 2008, 2009; Abbasi et al. 2011, 2012;
origin models invoke strangelet production in molecular
- Amenomori et al. 2012; Aartsen et al. 2014; Bartoli
o clouds (Kotera et al. 2013) or in neutron stars (Perez-
et al. 2013; Abeysekara et al. 2014; Bartoli et al. 2015).
r Garcia et al. 2014).
t For a recent summary of experimental results, we refer
s to Di Sciascio & Iuppa (2014). The statistics and reso- The measurement of cosmic-ray anisotropies with rel-
a ative intensity below 10−3 is an observational challenge,
lution of these experiments allow for a two-dimensional
[
since it is necessary to account for minuscule variations
reconstruction in the form of anisotropy sky maps. The
2 studies have revealed significant anisotropies at both in the acceptance and uptime of the detector carrying
v large and small angular scales. At the largest scales, out the measurement. For illustration, Fig. 1 shows a
7 the anisotropy is approximately dipolar and has a rela- simulated realization of the cosmic-ray anisotropy fol-
7 tive intensity on the order of 10−3. An explanation of lowing the model of Ahlers (2014). This model was cho-
8 senduetoitsqualitativeagreementwiththeobservation
theamplitudeandphaseofthedipoleanisotropyischal-
7 of anisotropy in TeV cosmic rays by IceCube (Aartsen
lenging, buttheobservationsarequalitativelyconsistent
0 et al. 2014) and HAWC (Abeysekara et al. 2014). The
with diffusive propagation of cosmic rays from Galactic
. cosmic-ray anisotropy is expected to remain constant in
1 sources(Ptuskinetal.2006;Erlykin&Wolfendale2006;
the celestial coordinate system over the period of the
0 Blasi & Amato 2012; Mertsch & Funk 2015).
observation. Ground-based cosmic observatories with a
6 Medium- and small-scale anisotropies have been ob-
1 served at the 10−4 level. These features are less under- limited field of view, however, are exposed to different
: stood but could be a combined effect of nearby cosmic- parts of the celestial sphere as the Earth rotates during
v one sidereal day. As an example, the sky map of Fig. 1
raysources(Salvati&Sacco2008;Biermannetal.2013)
i
X and the local interstellar magnetic field structure, which indicates the instantaneous field of view of the HAWC
detector (Abeysekara et al. 2014) (at latitude 19◦ N) at
r canintroduceanenergy-dependentmagneticmirrorleak-
a age (Drury & Aharonian 2008), preferred cosmic-ray alocalsiderealtimeof9h. Hence,theobservedeventdis-
tribution accumulated over many sidereal days depends
transport directions (Malkov et al. 2010), or magnetic
not only on the cosmic-ray anisotropy but also on the
lenses (Battaner et al. 2011, 2015). The observed power
nonuniform and time-dependent detector exposure. It is
spectrum of cosmic-ray anisotropies agrees well with the
therefore necessary to reconstruct a reference map that
expectedeffectofcosmic-rayscatteringinthelocalmag-
represents the response of the detector to an isotropic
netic field (Giacinti & Sigl 2012; Ahlers 2014; Ahlers &
cosmic-ray flux. For ground-based detectors, in which
1WIPAC & Department of Physics, University of Wisconsin– the atmosphere acts as part of the observatory, suffi-
Madison,Madison,WI53706,USA ciently accuratesimulations ofthe detectorexposureare
2DepartmentofPhysics&Astronomy,UniversityofRochester, usually not achievable. Therefore, most analyses use the
Rochester,NY14627,USA data themselves to estimate the relative intensity and
3CentroUniversitariodelosValles,UniversidaddeGuadalajara,
detector exposure simultaneously.
Guadalajara,Jalisco44130,M´exico
4nowatDepartmentofPhysics,UniversityofMaryland,College Examples of this technique are the time-
Park,MD,USA scrambling (Alexandreas et al. 1993) and direct-
†[e-mail] integration (Atkins et al. 2003) methods, in which the
2 Ahlers, BenZvi, Desiati, D´ıaz–V´elez, Fiorino & Westerhoff
anisotropy tion 5 before concluding in Section 6.
North
60◦ 2. MAXIMUM LIKELIHOOD METHOD
Inthefollowing,wewillassumethatthetotalaccumu-
Zenith 30◦ lated detector exposure E can be expressed as a product
East West of its angular-integrated exposure E and relative accep-
300◦ 240◦ 180◦ 120◦ 60◦ tanceAintermsofazimuthangleϕ(fromnorthincreas-
ing to the east) and zenith angle θ as
−30◦
E(t,ϕ,θ)(cid:39)E(t)A(ϕ,θ). (1)
South
−60◦ Equatorial Without loss of generality, we require that the relative
(cid:82)
acceptance is normalized to dΩA(ϕ,θ) = 1. This ap-
proximation assumes that the relative acceptance of the
-0.0008 0.0008
detectorremainsapproximatelyconstantovertime. The
Figure 1. Simulatedcosmicrayanisotropyinequatorialcoordi-
ansatzisidenticaltotheapproachusedindirectintegra-
natesusingthemodelofAhlers(2014). Forillustration,weindicate
the instantaneous field of view of the HAWC observatory (at lat- tion or time scrambling.
itude 19◦ north) at a local sidereal time of 9h and a zenith angle Let us also assume that the flux of cosmic rays at
cut at 60◦. The time-integrated field of view corresponds to the
the energies of interest remains constant as a function
declinationrange−41◦<δ<79◦.
of time, varying only as a function of celestial longitude
α (right ascension) and latitude δ (declination). Due to
rate of events observed in a detector as a function of the strong diffusion of cosmic rays in the Galactic envi-
local sidereal time is integrated against the relative ronment,thefluxisdominatedbyanisotropictermφiso.
acceptance of the detector during an integration period Hence, the total flux can be expressed as
(or scrambling interval) ∆t. The idea behind this
method is that variations in the event rates introduced φ(α,δ)=φisoI(α,δ), (2)
by the cosmic-ray anisotropies will average out to some
where I(α,δ) is the relative intensity of the flux as a
extent as the detector observes different parts of the
functionofpositioninthesky. Theanisotropyisdefined
celestial sphere over the course of each sidereal day. The
as the deviation δI = I −1 (cid:28) 1. Note that this ansatz
resultisanestimateofthenumberofeventsexpectedin
ignores anisotropies associated with the relative motion
the detector between time t and t+∆t. Subtracting the
oftheEarthwithrespecttotheSun. Wewillcomeback
expected number of events from the actual observations
to this subtlety in the discussion section.
yields a residual counts map which can be explored for
The local horizontal coordinate system and the ce-
anisotropy.
lestial (or equatorial) coordinate system are related
The integration time interval ∆t acts as an effec-
via a time-dependent transformation. We define n =
tive smoothing parameter for the counts map, since
(cosαcosδ,sinαcosδ,sinδ) as the unit vector corre-
the method averages cosmic-ray arrival directions over
sponding to the coordinates (α,δ) in the right-handed
angular scales of 15◦(∆t/1hr). In principle, choosing
equatorialsystem. Similarly,theunitvectorcorrespond-
∆t = 24h would produce a residual counts map with
ing to the coordinates (θ,ϕ) in the right-handed local
featurescoveringthefullsky(360◦). However,fordetec-
system is n(cid:48) = (cosϕsinθ,−sinϕsinθ,cosθ). The two
tors located in the middle latitudes, the instantaneous
unit vectors are related via a time-dependent coordinate
exposure of the detector does not match the full daily
transformation n(cid:48) = R(t)n. For an experiment located
exposure, since the 24h integrated field of view is much
atgeographiclatitudeΦandlongitudeΛ(measuredeast
largerthantheinstantaneousfieldofview,cf.Fig1. Asa
from Greenwich), the transformation is
result, large-scale structures in the residual counts map,
in particular the dipole, are strongly attenuated when (cid:32)−cosωtsinΦ −sinωtsinΦ cosΦ(cid:33)
using these methods. R(t)= sinωt −cosωt 0 , (3)
To improve estimates of large-scale anisotropy using cosωtcosΦ sinωtcosΦ sinΦ
detectorsinthemiddlelatitudes,wedescribeamaximum
likelihood construction that can be used to disentangle whereω =2π/24handthelocalsiderealtimetisrelated
the anisotropy from detector effects. The technique is to the sidereal time at Greenwich t(cid:48) by t=t(cid:48)+Λ/ω.
based on the same ansatz used by the time-scrambling To simplify calculations on the local and celestial
ordirect-integrationmethod,thatthetotalaccumulated spheres, the sky is binned into pixels of equal area
exposure of the detector can be factorized into a time- ∆Ω using the HEALPix parametrization of the unit
dependenteventrateandatime-independentrelativeac- sphere(Gorskietal.2005). Tomaketheequationsmore
ceptance map. We begin by describing the technique transparent, we use roman indices for pixels in the lo-
in Section 2. In Section 3, we apply the maximum- cal sky map and fraktur indices for pixels in the celestial
likelihoodmethodtosimulateddataandshowthatlarge- sky map. Time bins are indicated by greek indices. For
and small-scale anisotropies can be reconstructed with instance, the data observed at a fixed sidereal time bin
relatively little of the distortion observed in the direct- τ can be described in terms of the observation in local
integration or time-scrambling techniques. We compare horizontal sky with bin i as n or transformed into the
τi
our method to alternative techniques in Section 4. We celestial sky map with bin a as n .
τa
then discuss analysis methods of large- and small-scale Consider an angular element of the local coordinate
anisotropiesofthereconstructedanisotropymapsinSec- sphere ∆Ω corresponding to coordinates (θ ,ϕ ). The
i i i
Maximum-Likelihood Technique for Reconstructing Cosmic-Ray Anisotropy 3
number of cosmic rays expected from this location in a Because of this degeneracy, we must choose whether to
sidereal time interval ∆t with central value t is account for the local excess of cosmic rays as originating
τ τ
in an anisotropy in their relative intensity or originating
µ (cid:39)I N A , (4)
τi τi τ i in a variation in the local acceptance.
where N ≡ ∆t φisoE(t ) gives the expected number A natural choice is that the anisotropy is normalized
τ τ τ (cid:82)
of isotropic background events in sidereal time bin τ. to dαδI(α,δ)=0foralldeclinationsδ,consistentwith
(cid:82)
The quantity A ≡ ∆Ω A(θ ,ϕ ) is the binned relative the definition dΩδI(α,δ) = 0. This condition can also
i i i i
acceptance of the detector for angular element i, and be formulated in terms of a spherical harmonics expan-
I ≡ I(R(t )n(cid:48)(Ω )) is the relative intensity observed sion of the relative intensity in the equatorial coordinate
τi τ i
in the local horizontal system during time bin τ. For system (α,δ), as pointed out by Iuppa & Di Sciascio
simplicity, all expressions that follow assume equal bin (2013). In general, the relative intensity can be decom-
sizes∆Ω=4π/Npixonthelocalandcelestialspheresand posed as a sum over spherical harmonics Y(cid:96)m,
N equalsiderealbinswith∆t=24h/N . However,
time time
a more general binning is also possible. (cid:88) (cid:88)(cid:96)
Given µτi, the likelihood of observing n cosmic rays is δIa = (cid:98)a(cid:96)mYa(cid:96)m. (12)
given by the product of Poisson probabilities (cid:96)≥1m=−(cid:96)
L(n|I,N,A)=(cid:89)(µτi)nτie−µτi , (5) Othuercnonordmitaiolinzaation=co0nfdoirtiaolnl (cid:96)c.anTthhisenprboejeecxtipornesssiegdniafis-
n ! (cid:98)(cid:96)0
τi τi cantly reduces the reconstruction of the low-(cid:96) multipole
where n is the number of events observed in the local components of the anisotropy, as we will discuss in the
τi
pixel i during time bin τ. This likelihood can be max- following sections.
imized to provide estimators of the relative acceptance Note that the true multipole moments (cid:98)a(cid:96)m are an (in-
functionA andtheexpectedisotropicbackgroundcount finite) superposition of the pseudo multipole moments
i
Nτ. a(cid:96)m, which are defined as in Eq. (12), but for the prod-
Consider the case of the null hypothesis of no uct of the relative intensity with the weight function w
anisotropy, which we denote I(0) =1. Given the bound- of the field of view. Provided that the weight function
(cid:80) a is azimuthally symmetric, w(α,δ)=w(δ), the true mul-
ary condition A = 1, the maximum likelihood esti-
i i tipole moments a are a linear superposition of pseudo
mators of A and N are (cid:98)(cid:96)0
i τ multipole moments a . In practice, we can hence use
(cid:96)(cid:48)0
Nτ(0) =(cid:88)nτi, (6) athe=no0rmfoarliazlalt(cid:96)i.onInctoenrdmitsioonftah(cid:96)e0b=inn0edforrelaaltliv(cid:96)etionteennssuitrye
i (cid:98)(cid:96)0
A(0) =(cid:88)n (cid:46)(cid:88)n . (7) a(cid:80)ndwweYig(cid:96)h0δtIfu=nct0iofonrtahlils(cid:96).is equivalent to the condition
i τi κj a a a a
τ κj
2.2. Maximum Likelihood Algorithm
To allow for the possibility of anisotropy, we maximize
the likelihood ratio The maximum (I(cid:63),N(cid:63),A(cid:63)) of the likelihood ratio (8)
must obey the implicit relations
L(n|I,N,A)
λ= (8) (cid:88) (cid:46)(cid:88)
L(n|I(0),N(0),A(0)) I(cid:63) = n A(cid:63) N(cid:63), (13)
a τa κa κ
of signal over null hypothesis in N, A, and I. τ κ
(cid:88) (cid:46)(cid:88)
N(cid:63) = n A(cid:63)I(cid:63) , (14)
2.1. Invariance under Declination-Scaling τ τi j τj
i j
We will demonstrate in the following that the
(cid:88) (cid:46)(cid:88)
maximum-likelihood method can be used to search for A(cid:63) = n N(cid:63)I(cid:63) , (15)
i τi κ κi
anisotropyonallangularscales. However,beforewecon-
τ κ
tinue,notethateventsrecordedinafixedpositioninthe
local coordinate system can only probe the cosmic-ray together with (cid:80)awaYa(cid:96)0δIa(cid:63) = 0 and (cid:80)iA(cid:63)i = 1. In
flux in a fixed declination band δ. Hence, the expecta- Eq. (13), we introduced the binned quantity Aτa ≡
tion values (4) are invariant under the rescaling ∆ΩaA(RT(tτ)n(Ωa)) corresponding to the relative ac-
ceptanceseenintheequatorialcoordinatesysteminpixel
I → I(cid:48) ≡I/a(δ)/b, (9)
a during time bin τ.
N →N(cid:48)≡Nbc, (10) Equations (13) to (15) correspond to a nonlinear set
A→A(cid:48)≡Aa(δ)/c, (11) of equations that cannot be solved in an explicit form.
However, one can approach the best-fit solution via the
where a(δ) is an arbitrary function of declination and following iterative method:
the normalization factors b and c are defined such that
(cid:80) A (cid:48) = 1 and (cid:80) δI(cid:48) = 0 for the new values. In (i) Initialize at the maximum of the null hypothesis,
τ τ a a
other words, the maximum-likelihood method is sen- (I(0),N(0),A(0)).
sitive to anisotropy in right ascension but is insensi-
tive to variations in intensity across declination bands. (ii) Evaluate I(n+1) by inserting (I(n),N(n),A(n)) into
This is also known to be a limitation of other recon- the right-hand side of Eq. (13).
struction methods, like direct integration or time scram-
bling(Amenomorietal.2005;Iuppa&DiSciascio2013). (iii) Remove the m = 0 (pseudo) multipole moments
4 Ahlers, BenZvi, Desiati, D´ıaz–V´elez, Fiorino & Westerhoff
of δI(n+1), i.e., in the equatorial coordinate system 3
(cid:80) w Y(cid:96)0δI(n+1) →0.
a a a a 2
(iv) Evaluate N(n+1) by inserting (I(n+1),N(n),A(n)) 4]−
0
1 1
into the right-hand side of Eq. (14). [
1
−
(v) EvaluateA(n+1) byinserting(I(n+1),N(n+1),A(n)) ()0 0
into the right-hand side of Eq. (15). N
/ 1
()n −
(vi) Renormalize N(n+1) and A(n+1) as N(n+1) → N
2
N(n+1)c and A(n+1) → A(n+1)/c with normaliza- −
tion factor c=(cid:80) A(n+1). 3
i i − 0 3 6 9 12 15 18 21 24
localsiderealtime[h]
(vii) Repeatfromstep(ii)untilthesolutionhassufficient
300
convergence, i.e., the likelihood ratio in Eq. (8) has
∆χ2 (cid:39)2ln(λ(n+1)/λ(n))(cid:28)1.
250
Note that the cosmic-ray anisotropy obtained in the
first iteration step, δI(1), corresponds to the result that ()1)200
would be obtained by the method of direct integra- /λ
tion (Atkins et al. 2003). This method was also used ()n 150
λ
in the recent HAWC analysis (Abeysekara et al. 2014). (n
2l100
Thesuccessiveiterationstepsofthemaximum-likelihood
method reoptimize the relative acceptance A and the
50
isotropic background rate N for a given anisotropy. We
will study this optimization process more quantitatively
0
in the following section. 5 10 15 20
iterationstep
3. SIMULATION AND PERFORMANCE Figure 2. Toppanel: Therelativeoptimizationoftheisotropic
To demonstrate the reconstruction of a cosmic-ray expectationvalueintermsofN(n)/N(0)−1for20iterationsteps
(light to dark colors). Bottom panel: The progressive log-
anisotropy with the maximum likelihood technique,
likelihood values of the iteration. Note that the method already
we simulated a set of arrival directions based on the convergesafterabout10iterationsteps.
anisotropy shown in Fig. 1. The simulation is based on
a random realization of a relative intensity δI, which
follows a power-law spectrum with (cid:96) ≥ 1 of the form time bins). In each stable period, we assume that Nτ
C = 10−7(18/(2(cid:96)+1)/((cid:96)+1)/((cid:96)+2)), after the model has a normal fluctuation of 1% around the background
(cid:96)
ofAhlers(2014). Thismodelwaschosenduetoitsqual- expectation. This model mimics the data variation of
itative agreement with the observation of anisotropy in the analysis by Abeysekara et al. (2014).
TeV cosmic rays by IceCube (Aartsen et al. 2014) and Theresultsofthereconstruction6 areshowninFigs.2
HAWC (Abeysekara et al. 2014). The celestial sphere and3. ThetoppanelofFig.2showstheoptimizationof
isbinnedfollowingtheHEALPixparametrization(Gorski the reconstructed background expectation N(n) as the
et al. 2005) with parameter nside = 64 corresponding to relative quantity N(n)/N(0) − 1. The iteration renor-
Npix =49,152pixelswithabinsizeofabout1◦ diameter. malizes N(n) at the level of 10−4 depending on the local
The simulated detector is located at geographic coor-
siderealtime. Therelativefeaturesintroducedbytheop-
dinates (Φ = 19◦N and Λ = 97◦W), the location of the
timization process are easy to understand. For instance,
HAWC observatory (Abeysekara et al. 2014). The in-
the bump at 9h local sidereal time corresponds to the
stantaneous field of view is restricted to zenith angles
instantaneous field of view indicated in Fig. 1. In com-
below 60◦. A projection of the instantaneous field of
parison with the simulated anisotropy shown in Fig. 1,
view at 9h local sidereal time onto the equatorial coor-
one can notice that at this time the detector observes a
dinate system was already shown in Fig. 1. We assume
part of the sky with a strong underfluctuation. The ini-
that the relative detector acceptance at any time follows
tial estimate N(0) is therefore too low and the iteration
A(θ,ϕ)∝cosθ[1+Asinθsin2(ϕ−ϕ )]withA=0.2and
0 compensates for this effect. The bottom panel of Fig. 2
ϕ = 10◦. The local acceptance maps are reconstructed
0 shows the consecutive log-likelihood values of the itera-
with the same resolution as the anisotropy map.
tion. Note that for this simulated data set the iteration
The expected isotropic event number is binned in
only requires about 10 steps to converge. We run the
N = 360 sidereal time bins with 4 min bin size and
time reconstruction for 10 more iteration steps to verify the
normalized such that (cid:80) N (cid:39) 5 × 1010, correspond-
τ τ convergence and stability of the method.
ing to the integrated event number of the recent HAWC
The map in Fig. 3a shows the expected relative in-
analysis (Abeysekara et al. 2014). We introduce a sim- tensity smoothed with a top-hat function with 10◦ ra-
plestatisticaltoymodelforthevariationoftheexpected
dius after application of the m = 0 filter. The maps
background rate N ; the accumulated data are expected
τ
to be stable over periods sampled from an exponential 6 WeimplementedtheiterativealgorithmasaPythonprogram,
distribution with expectation value of 20 min (5 sidereal whichwecanprovidetotheinterestedreaderuponrequest.
Maximum-Likelihood Technique for Reconstructing Cosmic-Ray Anisotropy 5
(a)
anisotropywitha(cid:96)0→0and10◦smoothing
60◦
30◦
300◦ 240◦ 180◦ 120◦ 60◦
−30◦
−60◦
Equatorial
-0.0006 0.0006
(b) (c)=(b)−(a)
reconstructedanisotropy(10◦smoothed)/iteration1 residualanisotropy(10◦smoothed)/iteration1
60◦ 60◦
30◦ 30◦
300◦ 240◦ 180◦ 120◦ 60◦ 300◦ 240◦ 180◦ 120◦ 60◦
−30◦ −30◦
−60◦ −60◦
Equatorial Equatorial
-0.0006 0.0006 -0.0003 0.0003
(d) (e)=(d)−(a)
reconstructedanisotropy(10◦smoothed)/iteration20 residualanisotropy(10◦smoothed)/iteration20
60◦ 60◦
30◦ 30◦
300◦ 240◦ 180◦ 120◦ 60◦ 300◦ 240◦ 180◦ 120◦ 60◦
−30◦ −30◦
−60◦ −60◦
Equatorial Equatorial
-0.0006 0.0006 -0.0003 0.0003
Figure 3. Panel(a): Expectedrelativeintensitywithm=0multipolefilter(seetext)and10◦smoothing. Panel(b/d): Reconstructed
relative intensity for the first and last iteration step, respectively. Panel (c/e): Residual anisotropy after subtracting the expected
anisotropyinpanel(a)fromthereconstructedanisotropyinpanels(b)and(d),respectively.
in Figs. 3b and 3d show the reconstructed anisotropy in lated to the Poisson variation of the event rate. We can
the first iteration step (corresponding to the result from makethisstatementmorequantitativeviaapowerspec-
direct integration) and after 20 iterations, respectively. trum analysis of expected and reconstructed anisotropy
ComparingthemapsinFigs.3band3aonenoticesthat, maps. The relative intensity can be decomposed as a
qualitatively, the expected small-scale features are al- sum over spherical harmonics, as in Eq. (12). Unfortu-
ready reproduced in the first iteration step. However, nately, due to the limited integrated field of view, the
the difference map in Fig. 3c indicates, that the residual true coefficients a cannot be unambiguously recon-
(cid:98)(cid:96)m
map has large-scale features that are misreconstructed. structed. However, for the present discussion of residual
On the other hand, the corresponding difference map in anisotropies in the iterative method, it is sufficient to
Fig. 3e after 20 iterations is closer to the expected map study the pseudo multipole moments, a , correspond-
(cid:96)m
in Fig. 3a. ing to the harmonic expansion of the anisotropy multi-
The residual anisotropy map shown in Fig. 3e is re- plied by the weight function of the field of view. In our
6 Ahlers, BenZvi, Desiati, D´ıaz–V´elez, Fiorino & Westerhoff
10−7 mal relative detector acceptance (7) and the isotropic
iteration20 background expectation (6) of the null hypothesis (δI =
iteration1 0) are equivalent to the estimates via direct integration
expected with an integration period of ∆t = 24 h. The map of
Fig. 3b shows the resulting anisotropy estimate of this
m
ectru 10−8 meAtnhootdh.er two-dimensional reconstruction technique is
p
s
er the Forward-Backward method used by Milagro (Abdo
w
o et al. 2009), which is closely related to the one-
p
o dimensional East-West method (Bonino et al. 2011).
d
pseu 10−9 These methods analyze the relative right-ascension
derivative of event rates, ∂ n/n, at each sidereal time,
α
either for individual declination bands (in the case of
Forward-Backward)ortheentirefieldofview(inthecase
Nevents=5×1010&θzenith<60◦ ofEast-West). Theanisotropycanthenbereconstructed
10−10 2 4 6 8 10 from the first derivative, noting that ∂αδI (cid:39) ∂αn/n, up
multipolemoment(cid:96) to an overall normalization constant in each declination
band. This again reflects the invariance (9) to (11) and
Figure 4. Pseudo power spectra of reconstructed and expected
anisotropymaps. Theerrorbarsonthepowerofthereconstructed theinabilitytoreconstructthea(cid:96)0momentsinequatorial
anisotropy is indicating the variance from Poisson statistics esti- coordinates.
matedfrom100runswiththesameanisotropyanddetectormodel.
Ingeneral,theForward-Backward(andtheEast-West)
Notethatthesmallbiasofthereconstructedpseudopowertoward
highervaluesfor(cid:96)≥6isduetothenoiselevelofthemapsaswe method has the advantage that the quantity ∂αn/n
willshowinSection5. does not explicitly depend on the instantaneous back-
ground level. For instance, the local deficit of back-
example, the weight function is simply equal to 1 for de- ground events at a local sidereal time of 9h in the ex-
clinations −41◦ < δ < 79◦ and 0 otherwise. From these ample shown in Fig. 1 does not affect the estimate
pseudo multipole moments, we can compute the pseudo of the derivative ∂ n/n. However, the reconstruction
α
power spectrum, methods discussed in Abdo et al. (2009) assume that
the relative detector acceptance is quasi-symmetric un-
1 (cid:88)
C = |a |2. (16) der Forward-Backward (or East-West) reflection, corre-
(cid:96) 2(cid:96)+1 (cid:96)m sponding to the transformation (θ,ϕ) → (θ,−ϕ). This
m
condition can be expressed in terms of the asymmetry
Figure 4 shows the pseudo power spectrum of the
(cid:15)≡[A(θ,ϕ)−A(θ,−ϕ)]/[A(θ,ϕ)+A(θ,−ϕ)]as|(cid:15)|(cid:28)1.
anisotropy maps for the first and last iteration step in
Afterproperrenormalizationofthereconstructedderiva-
comparison to the expected spectrum corresponding to
tive terms (Abdo et al. 2009) the leading order effect is
the true anisotropy with a = 0 and multiplied by the
(cid:98)(cid:96)0 an O((cid:15)2) correction of the anisotropy amplitude. We
weight function. To estimate the variance of the re-
expect that the condition |(cid:15)| (cid:28) 1 is met by most ob-
constructed power spectra due to Poisson statistics, we
servatories. However, our method is also applicable for
repeat the analysis 100 times with low-resolution maps
large asymmetries approaching (cid:15)2 (cid:39)1.
(n =16 with bin size ∆θ (cid:39)4◦) and show the central
side Iterative methods for the reconstruction of the
68% range of the data. One can see that the first it-
anisotropy have also been developed for and applied to
eration step drastically underestimates the power of the
data of Tibet-ASγ (Amenomori et al. 2005, 2010, 2012)
dipole ((cid:96)=1) and quadrupole ((cid:96)=2). This was already
and ARGO-YBJ (Bartoli et al. 2015). The Equi-Zenith
qualitatively visible in the map of Fig. 3c. On the other
Angle method (Cui & Yan 2003) estimates the isotropic
hand, the pseudo power spectrum of the last iteration
cosmic-ray background at a celestial bin a to be the av-
step agrees well with the low-(cid:96) multipoles. One can also
erage of those acceptance-corrected events that arrived
noticeasmallbiasofthe(cid:96)≥6multipolestowardslarger
from the same zenith angle band as the field of view
values. This is due to the noise level of the maps as we
wraps around the celestial equator. This corresponds
will discuss in detail in Section 5. to a generalization of the ansatz (1) to E = (cid:80) EsAs,
s
where the sum runs over the different equi-zenith-angle
4. COMPARISON TO OTHER METHODS
sectors with individual background rates Es and rela-
Various reconstruction methods have been developed tive acceptance As. It is straightforward to generalize
previously to extract the anisotropy from the isotropic our likelihood-based method with this ansatz (see Ap-
background. We have already mentioned the meth- pendix A). The likelihood-based iteration method pre-
ods of direct integration Atkins et al. (2003) and time sented in this paper has the advantage that it can be
scrambling (Alexandreas et al. 1993) in the introduc- derived from a firm statistical approach. If the relative
tion,whichhavebeenusedincosmic-rayanisotropystud- detector acceptance can be regarded as stable in local
ies of Super-Kamiokande (Guillian et al. 2007), Tibet- sidereal time over the entire field of view, it provides an
ASγ (Amenomori et al. 2006), Milagro (Abdo et al. even simpler iterative reconstruction method.
2008), IceCube (Abbasi et al. 2011, 2012; Aartsen et al.
5. HARMONIC ANALYSIS
2014), IceTop (Aartsen et al. 2013), HAWC (Abey-
sekaraetal.2014),andARGO-YBJ(Bartolietal.2013). We now turn to the harmonic analysis of the
The time-integration method is closely related to the anisotropy. We have already introduced the harmonic
maximum-likelihood method presented here. The opti- expansion of the relative intensity in Eq. (12). Of par-
Maximum-Likelihood Technique for Reconstructing Cosmic-Ray Anisotropy 7
ticular relevance for the theory of cosmic-ray diffusion is (a)
the strength of the dipole components ((cid:96)=1), but note large-scaleanisotropy(truncated,(cid:96) 3)/iteration20
≤
thatTeV–PeVcosmic-raydataalsoshowsignificantmul-
tipolemomentsatsmallerangularscales,likequadrupole 60◦
((cid:96)=2), octupole ((cid:96)=3), etc.
We have already shown in Section 3 that traditional 30◦
anisotropy methods can significantly underestimate the
low-(cid:96) pseudo power spectrum of the data. This effect 300◦ 240◦ 180◦ 120◦ 60◦
can be compensated by the iterative method presented
in Section 2. However, this method will not compen- −30◦
sate for the loss in power of the m = 0 coefficients of
the analysis. The observed spectrum is hence always a −60◦
systematic underestimation of the true anisotropy.
An additional uncertainty comes from the limited in-
tegrated field of view of most observatories. In the ideal -0.0005 0.0005
case of a 4π sky coverage, the multipole moments a of (b)
(cid:96)m
thereconstructedanisotropywouldcarryalltheinforma- small-scaleanisotropy((cid:96)≤3removed,10◦smoothed)/iteration20
tion of the anisotropy (except a ). However, as already
(cid:98)(cid:96)0
discussed earlier, the partial sky coverage of individual 60◦
experiments does only allow reconstructing the pseudo
multipolemomentsspectrumofthe(reduced)anisotropy 30◦
multiplied by the weight function.
The pseudo multipole moments a are related to the
(cid:96)m 300◦ 240◦ 180◦ 120◦ 60◦
true multipole moments a via a linear transformation
(cid:98)(cid:96)m
(see, e.g., the review by Efstathiou (2004))
−30◦
(cid:88)
a = K a , (17)
(cid:96)m (cid:96)m(cid:96)(cid:48)m(cid:48)(cid:98)(cid:96)(cid:48)m(cid:48) −60◦
(cid:96)(cid:48)m(cid:48) Equatorial
where the coupling matrix K depends on the multipole
spectrum b of the weight function of the field of view -0.0004 0.0004
(cid:96)m
(see Appendix B). For an azimuthally symmetric weight (c)
function the strength of the mixing between moments significance(unitsofσ,10◦smoothed)/iteration20
a and a is determined by the moments b with
(cid:96)m (cid:98)(cid:96)(cid:48)m k0
|(cid:96)−(cid:96)(cid:48)|≤k ≤(cid:96)+(cid:96)(cid:48). Thismixingcanbesmallforindivid- 60◦
ual moments as pointed out by Denton & Weiler (2015).
However, it is important to emphasize that the full mul- 30◦
tipole spectrum cannot be unambiguously reconstructed
fromapartialskycoveragesincetheinfinite-dimensional 300◦ 240◦ 180◦ 120◦ 60◦
matrix K is not invertible.
−30◦
5.1. Large-Scale Anisotropy
Whereas the full transition matrix K cannot be in- −60◦
Equatorial
verted, we can attempt an approximate reconstruction
of the low-(cid:96) anisotropy via a truncation of the multipole
-9 9
expansion after a maximum (cid:96). The corresponding trun-
cated matrix K(cid:48) is then invertible. For instance, assum- Figure 5. Panel (a): Reconstructedlarge-scaleanisotropyfrom
thesolutionofthematrixequationEq.(17)truncatedafter(cid:96)=3.
ing a pure dipole anisotropy, (cid:96) ≤ 1, and a uniform sky
Panel (b): Residual small-scale anisotropy after subtracting the
coverage between declination δ1 and δ2 gives the dipole map in panel (a) from the full reconstructed anisotropy Panel
transition elements (c): Significance of the small-scale structure using Eq. (24). To
distinguish excesses from deficits we multiply by the sign of the
K(cid:48) = 1(cid:0)sin3δ −sin3δ (cid:1), (18) smoothedresidualanisotropyδIsmall.
1010 2 2 1
K1(cid:48)111 = 14(cid:0)3(sinδ2−sinδ1)+sin3δ1−sin3δ2(cid:1), (19) 7m9o◦mgeinvtinsgofKth1(cid:48)e11p1ro(cid:39)je0c.t9e2d.diHpeonlec,ea, thaenpdsaeudo,hmavuelttipooblee
11 1-1
and K(cid:48) =K(cid:48) . correctedbyamoderatefactor1/K1(cid:48)111 (cid:39)1.09torecover
1-11-1 1111 thetruemomentsaˆ andaˆ . Ontheotherhand,using
ThedipolestrengthperpendiculartotheEarth’srota- 11 1-1
the same zenith cut at the location of IceCube/IceTop
tion axis can then be estimated as
gives δ = −90◦ and δ = −30◦, yielding K(cid:48) (cid:39) 0.16
(cid:114) 3 (cid:112)|a |2+|a |2 and a c1orrection factor21/K(cid:48) (cid:39)6.4. 1111
A(cid:98)⊥ = 4π 11K(cid:48) 1-1 . (20) However,weemphasizeth1a1t1t1histreatmentisonlycor-
1111 rect under the assumption that the true anisotropy is
For instance, in our MC simulation with a zenith cut of dominated by a dipole. In Fig. 5a, we show a map of
60◦ andalatitudeΦ(cid:39)19◦,wehaveδ =−41◦ andδ = the reconstructed large-scale anisotropy including dipole
1 2
8 Ahlers, BenZvi, Desiati, D´ıaz–V´elez, Fiorino & Westerhoff
((cid:96) = 1), quadrupole ((cid:96) = 2), and octupole ((cid:96) = 3) with 6
a truncation of the matrix in Eq. (17) after (cid:96)=3. Note 1
−
tahlsaotatshseummeosdesilgonfifiAchalnetrsp(o2w0e1r4a)tuhseigdhfeorrmthuilstispiomleuslawtiiothn exp/C(cid:96) 0
(cid:96) > 3. The truncation after (cid:96) = 3 is hence not per se C(cid:96)
6
justified by this model. However, we can still use the 10−−6
reconstructed large-scale anisotropy map in Fig. 5a as a
backgroundmodeltodefinethesmall-scaleanisotropyin Nevents=5×1010&θzenith<60◦ expectedC(cid:96)
noiseN
the full anisotropy map via subtraction.
10−7 sampleC(cid:96)
pseudoC(cid:96)
5.2. Small-Scale Anisotropy
estimatedC(cid:96)
tuTrehseinsttahteisfitincaallrseigconnifisctraunccteedofmraespidisuaulsuaanlilsyotersotpimyafteead- ctrum 10−8 (cid:98)
e
p
using the method introduced by Li & Ma (1983) for ap- s
er
plicationsingamma-rayastronomy. Adirectapplication ow 10−9
of this method does not account for the optimization p
process of the time-dependent exposure. However, it is
rather straightforward to generalize the method of Li & 10−10
Ma (1983) to our case.
We begin by dividing the reconstructed relative inten-
ssictayleinfteoataurceosn,tIri=buItiloarngeo+f lIasrmgeal-ls.caFloerfeeaactuhrpesixaenldasinmathlle- 10−11 5 10 15 20 25 30
multipolemoment(cid:96)
celestialsky,wedefineexpectedon-sourceandoff-source
event counts in a disc of radius ψ centered on that pixel. Figure 6. Multipole power spectra of the final iteration map of
Given the set of pixels D , the observed and expected theexample. Thedashedlineshowsthemodelexpectationofthe
a inputpowerandtheopencirclesshowthesampledpowerspectrum
counts are
oftheinputanisotropymap. Theopenboxesarethepseudopower
(cid:88) (cid:88) spectrum for the reconstructed anisotropy. The filled diamonds
n = n , (21)
a τb showtheestimatedpowerspectrumaccountingforthenoiselevel
b∈Da τ (dottedline)andthefieldofviewusingthemethodsdescribedin
(cid:88) (cid:88) the main text. The upper plot shows the relative scatter of the
µ = A N I , (22) reconstructedpowerspectrumfromtheexpectedinputpower.
a,on τb τ b
b∈Da τ
(cid:88) (cid:88)
µ = A N Ilarge. (23) spectrum. However, in certain situations one can make
a,off τb τ b
additional assumptions about the ensemble-averaged ex-
b∈Da τ pectationvaluesofthemultipolecomponents. Inthefol-
The significance map (in units of Gaussian σ) is then lowing,wewillassumethattheharmoniccoefficientsare
calculated as Gaussian random fields, which in the ensemble-average
√ (cid:18) µ (cid:19)1/2 follow(cid:104)(cid:98)a(cid:96)m(cid:98)a∗op(cid:105)=δmpδ(cid:96)o(cid:104)C(cid:98)(cid:96)(cid:105). Inthisparticularcase,we
Sa = 2 −µa,on+µa,off +nalogµa,on . (24) can recover the ensemble-averaged power spectrum (cid:104)C(cid:98)(cid:96)(cid:105)
a,off via the relation
In Fig. 5b, we show the residual anisotropy map after (cid:88)
subtraction of the large-scale anisotropy map shown in (cid:104)C(cid:96)(cid:105)= M(cid:96)(cid:96)(cid:48)(cid:104)C(cid:98)(cid:96)(cid:48)(cid:105)+N(cid:96). (25)
Fig. 5a. The corresponding significance map using (24) (cid:96)(cid:48)
is shown in Fig. 5c. Usually, the significance is multi- The transfer matrix M is known from the study of
plied by the sign of the (smoothed) anisotropy (middle temperature anisotropies in the cosmic microwave back-
panel) to distinguish deficits (blue) from excesses (red), ground (Efstathiou 2004). However, for our situation of
andwefollowherethesameconvention. Wecanseethat cosmic-ray anisotropies we have again to account for the
after the subtraction of large-scale features there is still fact that the m = 0 moments are filtered out by the re-
significant power in small-scale features at a significance construction. This leads to a modified expression for M
level of 8σ. that we discuss and provide in Appendix B.
Another approach to study small-scale structure is via In Eq. (25), we have also introduced the noise power
the power spectrum (16), which quantifies the absolute spectrumN ,whichcaningeneralbecalculatedfromthe
(cid:96)
amplitudeofthemultipolecomponentsbutignorestheir relative intensity variance from the likelihood function
phases. We have already used the power spectrum to (5). In the following, we will use an approximation that
quantify the performance of our iterative method de- onlydependsonpixel-by-pixelPoissonnoise,whichgives
scribed in Section 2. We can also use this quantity in a flat spectrum N =N with
(cid:96)
situations where the small-scale anisotropy structure is
dominated by random scattering of cosmic rays in local N (cid:39) 1 (cid:88) wa2∆Ω2 . (26)
magneticfields(Giacinti&Sigl2012;Ahlers2014;Ahlers 4π (cid:80) n
& Mertsch 2015; L´opez-Barquero et al. 2015). a τ τa
As mentioned in the previous section, the true multi- In Fig. 6, we show the estimated power spectrum in-
pole spectrum cannot be recovered unambiguously from ferred from the relative intensity map of the last iter-
apartialskycoverage,andthesameistrueforthepower ation step shown in Fig. 3d. The dashed line corre-
Maximum-Likelihood Technique for Reconstructing Cosmic-Ray Anisotropy 9
sponds to the expected input power from which the in- tions in the primary cosmic-ray flux during the observa-
put power spectrum (black data points) is sampled as tionperiod. Thesolarpotential,withavariationontime
a Gaussian random field. The sampling introduces a scales of 11 years, only affects low-energy data and can
scatter of (∆C )2 = 2C2/(2(cid:96)+1) around the expected be neglected for cosmic rays in the TeV to PeV energy
(cid:96) (cid:96)
input power (cosmic variance). The green data points range(Munakataetal.2010). Ontheotherhand,therel-
showthepseudopowerspectrum, notaccountingforthe ativemotionoftheEarthinthesolarsystemisexpected
weight function. The horizontal blue dotted line shows tobevisibleasasolardipole(Compton&Getting1935).
the noise level of Eq. (26). Note that the noise level Thiseffectwouldbevisibleifthedataisbinnedinterms
scales as N ∝ 1/N with the total number of events. of solar time. Indeed, a solar dipole anisotropy has been
tot
The red data points are the best estimators of the true observed at the level of 10−4 in multi-TeV cosmic-ray
power spectrum (cid:104)C(cid:98)(cid:96)(cid:105). The variance of the pseudo and data (Amenomori et al. 2004, 2006; Abdo et al. 2009;
estimated power spectra are given in Appendix B. Abbasi et al. 2011, 2012; Bartoli et al. 2015). However,
the solar dipole is expected to average out for the data
6. SUMMARY binned in sidereal time as long as the observation period
covers an integer number of full years.
In this paper, we have discussed a novel two-
The method presented in this paper is designed for
dimensional cosmic-ray reconstruction method. It is
midlatitude observatories that are exposed to different
based on a maximum likelihood analysis that provides
parts of the celestial sky as the Earth rotates. The
implicitbest-fitexpressionsfortherelativeintensity,rel-
method is also applicable to the case of IceCube/IceTop
ativedetectoracceptance,andbackgroundexpectations.
located at the South Pole. However, for these observa-
Wehaveprovidedadetailediterativemethodonhowthe
tories, the instantaneous field of view is identical to the
relative intensity can be reconstructed from these im-
time-integrated one. As a result, the iteration does not
plicit maximum likelihood solutions. The performance
improve the estimates from direct integration or time
of this likelihood-based method was studied via a simu-
scrambling. On the other hand, the special location in
lated example, mimicking the position and performance
combination with a limited field of view makes these
of the HAWC observatory.
observatories particularly insensitive to low-(cid:96) multipoles
In general, ground-based observatories are insensitive
due to projection effects of the anisotropy discussed in
to cosmic-ray anisotropy variations across declination
the text.
bands. Intermsofasphericalharmonicexpansionofthe
In this paper, we have only considered the case of the
anisotropy in equatorial coordinates this corresponds to
dataanalysisofindividualobservatories. Agreatadvan-
a filter of m=0 multipole moments, which introduces a
tageofthelikelihood-basedreconstructionmethodisthe
systematic underestimation of the observed anisotropy.
straightforward generalization to combined anisotropy
In particular, the dipole anisotropy, which is a crucial
studiesofdatasetsfrommultipleobservatorieswithover-
observable for the study of cosmic ray diffusion in our
lapping integrated field of view, see e.g. (D´ıaz–V´elez, J.
Galaxy, can only be observed as a projection onto the
C. 2015). The expectation value (4) can then be sim-
celestial equator.
ply generalized to a sum over data sets with individual
Thishasimportantconsequencesfortheinterpretation
detector exposures but same anisotropy, as long as the
of experimental dipole data. If the dipole orientation
rigidity distributions of the data sets are very similar.
changes with rigidity the projected dipole can exhibit
We would also like to emphasize that cosmic-ray ob-
rigiditymodulationsintroducedbytheprojectionontop
servations via satellites like Fermi could have an advan-
of a true rigidity dependence of the dipole amplitude.
tage since the observatory can be tilted during the ob-
If the dipole vector aligns with the celestial poles these
servation. In this case, it is possible to break the de-
projection effects can lead to a drastic reduction of the
generacy between local acceptance effects and cosmic-
observed dipole accompanied by a phase-flip.
ray anisotropy. In principle, it is then possible to have
In addition, the limited integrated field of view of ob-
a full reconstruction of the anisotropy without projec-
servatories affects the reconstruction of the multipole
tion effects, provided that additional systematic effects
moments. As as consequence, even the m = 0 filtered
are under control.
anisotropycannotbeunambiguouslyreconstructed. This
can introduce a large systematic uncertainty, in partic-
ular, for the low-(cid:96) multipole moments of the anisotropy. Acknowledgements. The authors acknowledge sup-
We have discussed strategies how to account for this ef- portbytheNationalScienceFoundation(PHY-1306465,
fect in the multipole reconstruction. PHY-1306958,andPHY-1308033),bytheDepartmentof
We would like to conclude with a few remarks. In our Energy (DE-SC0008475) and by the Wisconsin Alumni
anisotropy ansatz (2), we have not accounted for varia- Research Foundation.
APPENDIX
A. GENERALIZATION OF THE EXPOSURE ANSATZ
We can generalize the ansatz (1) by expressing the total accumulated exposure E as a sum over disjoint sky sectors,
whose union covers the entire field of view. As before, we assume that the exposure in each sector can be expressed as
a product of its angular-integrated exposure Es and relative acceptance in terms of azimuth ϕ and zenith angle θ as
(cid:88)
E(t,ϕ,θ)(cid:39) Es(t)As(ϕ,θ). (A1)
sectors
10 Ahlers, BenZvi, Desiati, D´ıaz–V´elez, Fiorino & Westerhoff
The partition of the field of view into multiple sectors is at this stage arbitrary and should in general be guided by
the property of the data. For instance, the Equal-Zenith-Angle Method divides the sky into ring-like sectors with
limited zenith angle range, i.e., we have sector weight functions defined by ws(θ) = 1 if θ < θ < θ and otherwise
s s+1
ws(θ)=0. With this new ansatz, the best-fit relative acceptance and background rate of the null hypothesis become
Ns(0) =(cid:88)wsn , As(0) =(cid:88)wsn (cid:46)(cid:88)wsn . (A2)
τ i τi i i τi j κj
i τ κj
Here, we have introduced the weight function ws of the sector s which is equal to 1 if the pixel i is located in the
i
sector and 0 otherwise. The maximum of the signal hypothesis now obeys the implicit relation
(cid:88) (cid:46)(cid:88) (cid:88) (cid:46)(cid:88) (cid:88) (cid:46)(cid:88)
I(cid:63) = n As(cid:63)Ns(cid:63), Ns(cid:63) = wsn As(cid:63)I(cid:63) , As(cid:63) = wsn Ns(cid:63)I(cid:63) . (A3)
a τa κa κ τ i τi j τj i i τi κ κi
τ sκ i j τ κ
B. POWER SPECTRUM ESTIMATOR AND VARIANCE
Inthefollowing,wewillassumethattheweightfunctionisazimuthallysymmetric. Inthiscase,thetransferfunction
is block-diagonal K =δ Tm(w) with block elements defined via a sum over Wigner-3j coefficients,
(cid:96)m(cid:96)(cid:48)m(cid:48) mm(cid:48) (cid:96)(cid:96)(cid:48)
(cid:96)(cid:88)+(cid:96)(cid:48) (cid:114)(2(cid:96)+1)(2(cid:96)(cid:48)+1)(2k+1)(cid:18)(cid:96) (cid:96)(cid:48) k(cid:19)(cid:18)(cid:96) (cid:96)(cid:48) k(cid:19)
Tm(b)=(−1)m b . (B1)
(cid:96)(cid:96)(cid:48) k0 4π 0 0 0 m −m 0
k=|(cid:96)−(cid:96)(cid:48)|
For the ensemble-averaged multipole moments, we can evaluate the transfer matrix to
2(cid:96)(cid:48)+1(cid:88) (cid:18)(cid:96) (cid:96)(cid:48) k(cid:19)2 [T0 (w)]2
M = (2k+1)W − (cid:96)(cid:96)(cid:48) . (B2)
(cid:96)(cid:96)(cid:48) 4π k 0 0 0 2(cid:96)+1
k
Notethattheunfamiliarlastterminthepreviousequationaccountsfortheprojectionofthepseudoangularmomentum
onto m(cid:54)=0 terms.
The variance can be expressed as (cid:104)∆C ∆C (cid:105)=V1 +V2 +V3 with
(cid:96) (cid:96)(cid:48) (cid:96)(cid:96)(cid:48) (cid:96)(cid:96)(cid:48) (cid:96)(cid:96)(cid:48)
2 (cid:88) (cid:88)(cid:88)
V(cid:96)1(cid:96)(cid:48) = (2(cid:96)+1)(2(cid:96)(cid:48)+1) (cid:104)C(cid:98)k(cid:105)(cid:104)C(cid:98)k(cid:48)(cid:105)Tkm(cid:96)(w)T(cid:96)m(cid:48)k(w)T(cid:96)mk(cid:48)(w)Tkm(cid:48)(cid:96)(cid:48)(w), (B3)
m(cid:54)=0 k k(cid:48)
2 (cid:88) (cid:88)
V(cid:96)2(cid:96)(cid:48) = (2(cid:96)+1)(2(cid:96)(cid:48)+1) (cid:104)C(cid:98)k(cid:105)[T(cid:96)m(cid:96)(cid:48)(u)Tkm(cid:96)(w)T(cid:96)m(cid:48)k(w)+T(cid:96)m(cid:48)(cid:96)(u)Tkm(cid:96)(cid:48)(w)T(cid:96)mk(w)] , (B4)
m(cid:54)=0 k
1 (cid:88) (cid:18)(cid:96) (cid:96)(cid:48) k(cid:19)2
V3 = (2k+1)U , (B5)
(cid:96)(cid:96)(cid:48) 2π k 0 0 0
k
where u is the multipole coefficient of the distribution ∆Ωw2/(cid:80) n and U the corresponding power spectrum.
(cid:96)m a τ τa (cid:96)
Since the variance matrix of C is the same as for C −N, we can express the variance of the true spectrum as
(cid:96) (cid:96)
(cid:104)∆C(cid:98)(cid:96)∆C(cid:98)(cid:96)(cid:48)(cid:105)=M(cid:96)−k1M(cid:96)−(cid:48)k1(cid:48)(cid:104)∆Ck∆Ck(cid:48)(cid:105). (B6)
REFERENCES
Aartsen,M.G.,etal.2013,Astrophys.J.,765,55
—.2014,Braz.J.Phys.,44,415
Abbasi,R.,etal.2011,Astrophys.J.,740,16
—.2012,Astrophys.J.,746,33
Abdo,A.,Allen,B.,Aune,T.,etal.2008,Phys.Rev.Lett.,101,221101
—.2009,Astrophys.J.,698,2121
Abeysekara,A.U.,etal.2014,Astrophys.J.,796,108
Ahlers,M.2014,Phys.Rev.Lett.,112,021101
Ahlers,M.,&Mertsch,P.2015,Astrophys.J.Lett.,815,L2
Alexandreas,D.E.,Berley,D.,Biller,S.,etal.1993,Nucl.Instrum.Meth.A,328,570
Amenomori,M.,etal.2004,Phys.Rev.Lett.,93,061101
—.2005,Astrophys.J.,633,1005
—.2006,Science,314,439
—.2010,Astrophys.J.,711,119
—.2012,Astropart.Phys.,36,237
Atkins,R.W.,etal.2003,Astrophys.J.,595,803
Bartoli,B.,etal.2013,Phys.Rev.,D88,082001
—.2015,Astrophys.J.,809,90
Battaner,E.,Castellano,J.,&Masip,M.2011,A&A,527,A79
—.2015,Astrophys.J.,799,157
Biermann,P.,Becker,J.,Seo,E.-S.,&Mandelartz,M.2013,Astrophys.J.,768,124
Blasi,P.,&Amato,E.2012,JCAP,1201,011