A Bayesian analysis of the 69 highest energy cosmic rays detected by the Pierre Auger Observatory

1.8 MB·

by  Alexander Khanin

Mon.Not.R.Astron.Soc.000,1–14(2016) Printed12thJanuary2016 (MNLATEXstylefilev2.2)
A Bayesian analysis of the 69 highest energy cosmic rays
detected by the Pierre Auger Observatory
Alexander Khanin1(cid:63) and Daniel J. Mortlock1,2
6 1Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, U.K.
1 2Department of Mathematics, Imperial College London, London SW7 2AZ, U.K.
0
2
Accepted2016????????.Received2016????????;inoriginalform2016??????????
n
a
J
ABSTRACT
1
Theoriginsofultra-highenergycosmicrays(UHECRs)remainanopenquestion.Sev-
1
eralattemptshavebeenmadetocross-correlatethearrivaldirectionsoftheUHECRs

] with catalogs of potential sources, but no definite conclusion has been reached. We
E reportaBayesiananalysisofthe69eventsfromthePierreAugerObservatory(PAO),
H thataimstodeterminethefractionoftheUHECRsthatoriginatefromknownAGNsin
. theVeron-Cety&Veron(VCV)catalog,aswellasAGNsdetectedwiththeSwiftBurst
h
AlertTelescope(Swift-BAT),galaxiesfromthe2MASSRedshiftSurvey(2MRS),and
p
an additional volume-limited sample of 17 nearby AGNs. The study makes use of a
-
o multi-level Bayesian model of UHECR injection, propagation and detection. We find
r that for reasonable ranges of prior parameters, the Bayes factors disfavour a purely
t
s isotropic model. For fiducial values of the model parameters, we report 68% cred-
a ible intervals for the fraction of source originating UHECRs of 0.09+0.05, 0.25+0.09,
[ −0.04 −0.08
0.24+0.12, and 0.08+0.04 for the VCV, Swift-BAT and 2MRS catalogs, and the sample
−0.10 −0.03
1 of 17 AGNs, respectively.
v
5 Key words: cosmic rays – methods: statistical
0
3
2
0
. 1 INTRODUCTION withtheirsources.WhileUHECRsarechargedparticlesand
1
therefore experience magnetic deflection as they propagate,
0 Cosmic rays (CRs) are highly accelerated protons and
theyaresufficientlyenergeticthatthetotaldeflectionisex-
6 atomic nuclei, some of which enter the Solar system and
pectedtobe∼2to∼10deg(e.g.MedinaTancoetal.1998;
1 reach the Earth. They are the most energetic particles ob-
Sigletal.2004;Dolagetal.2005),sothatsomeinformation
: servedinnature,withenergiesintherange109eVto1021eV
v about their points of origin should be retained.
i (see e.g. Kotera & Olinto 2011, Letessier-Selvon & Stanev
X 2011 for reviews). Association of UHECRs with catalogs of potential
sourcesismadepossiblebythefactthatUHECRswithen-
r A number of open scientific issues remain with re-
a ergiesofE (cid:38)5×1019eV areexpectedtohavecomefroma
spect to CRs, in particular ultra-high energy cosmic rays
(UHECRs) with arrival energies E (cid:38)1019eV. The study limitedradiusof∼100Mpc.Thisradiusissometimescalled
arr
theGreisen-Zatsepin-Kuzmin(GZK)horizon,andarisesdue
ofUHECRsiscomplicatedbythefactthattheyexperience
anabruptcutoffintheirenergyspectrumat∼4×1019eV,so to the fact that UHECRs at those energies scatter off the
cosmicmicrowavebackground(CMB)radiationinaprocess
that only small samples areavailable. Thelargest currently
available sample is the 69 events with E (cid:62)5.5×1019eV knownastheGZKeffect(Greisen1966,Zatsepin&Kuzmin
arr
1966).ThemeanfreepathoftheGZKeffectathighenergies
recorded by the Pierre Auger Observatory (PAO) between
isafewMpcandtheenergylossineachcollisionis20-50%.
2004 January 1 and 2009 December 31 (Abreu et al. 2010).
The resultant attenuation is very rapid, and is the cause of
One open issue in the study of UHECRs is the ques-
thecutoffintheUHECRenergyspectrumobservedbyboth
tionoftheirsources.Anumberofcandidates,suchasactive
HiRes(Abbasietal.2008)andPAO(Abrahametal.2008).
galacticnuclei(AGNs)andpulsarshavebeenproposed,but
studies have not been conclusive (see e.g. Kalmykov et al. A number of attempts have been made to find correl-
2013 for a review). The question of UHECR origins can ations between UHECR arrival directions and catalogs of
be studied by attempting to associate the arrival directions possible sources. Cross-correlation studies have been con-
ducted with galaxy catalogs, such as the Two Micron All-
Sky Survey (2MASS) Redshift Survey (2MRS) (Abraham
(cid:63) E-mail:[email protected] et al. 2009; Abbasi et al. 2010), as well as specific types
(cid:13)c 2016RAS
2 A. Khanin & D. J. Mortlock
of objects such as active galactic nuclei (AGNs) (Abraham isonareexploredinAppendixB.WeuseaHubbleconstant
et al. 2007; Abraham et al. 2008; George et al. 2008; Pe’Er of H = 70km/s/Mpc where required to convert between
0
et al. 2009; Watson et al. 2011) and BL Lacertae objects redshifts and distances.
(BL LAcs) (Tinyakov & Tkachev 2001). Overall, no clear
consensushasbeenreached.Differentstudieshavereported
differentdegreesofcorrelation,dependingonthestatistical
2 DATA
approach,theUHECRsample,andthepopulationofsource
candidates that was used. The most significant correlation 2.1 UHECR sample
wasreportedbythePierreAugerCollaboration,betweenar-
rivaldirectionsofUHECRswithenergiesE (cid:62)5.7×1019eV The sample of UHECR events that was used in this ana-
lysis were the 69 highest energy events recorded at the
and the positions of nearby AGNs (Abraham et al. 2007).
PAO between January 2004 and November 2009, as doc-
The result was supported by Yakutsk data (Ivanov 2009),
umented in Abreu et al. (2010). These are the events with
butnotbyHiRes(Abbasietal.2008)ortheTelescopeArray
observed energies E above the threshold E (cid:62)E =
(Abu-Zayyad et al. 2012). A more recent analysis of a lar- obs obs thres
5.7×1019eV.
gerPAOsamplehasshownaweakercorrelationthanbefore
The PAO is a CR observatory located in Argentina, at
(Abreu et al. 2010).
a longitude of 69.5◦ W and a latitude 35.2◦ S. PAO is a
The lack of consensus on these issues is partly due to
hybrid observatory, which means that it uses both surface
thedifficultyofanalyzingsuchsmallsamplesizes.Giventhe
detection (SD) and fluorescent telescope detection (FD) of
smallsizeoftheUHECRdatasets,itisimportanttoutilize
UHECRs. The observatory has SD plastic scintillators of a
asmuchoftheavailableinformationaspossible.Thiscanbe
total area of 3000 km2 and 4 FD telescopes.
achievedbyadoptingaBayesianmethodology,thatinvolves
The PAO’s total exposure of this data-set is (cid:15) =
models of the relevant physical processes. The first steps tot
20,370km2sryr and its relative exposure per unit solid
to such a comprehensive Bayesian work have been made in
angle, d(cid:15)/dΩ, is illustrated in Figure 1. The relative expos-
the recent work of Watson et al. (2011) and Soiaporn et al.
ureisdirectlyproportionaltoPr(det|r),theprobabilitythat
(2013).
aUHECRwillbedetectedifitarrivesfromdirectionr,but
Watson et al. (2011) analysed the 27 events that were
(cid:82)
is normalized so that (d(cid:15)/dΩ)dΩ=(cid:15) .
analysed in Abraham et al. (2007), and derived a posterior tot
PAO measures UHECR arrival directions with an un-
for the fraction that originated from AGNs in the Veron-
certainty of ∼1deg and arrival energies with a relative un-
Cetty&Veron(VCV)catalog(V´eron-Cetty&V´eron2006).
certainty of ∼12% (Letessier-Selvon et al. 2014).
To do so, they used a two-component parametric model
characterizedbyasourcerateΓandabackgroundUHECR
rate R. The model assumed that the UHECR arrival direc-
2.2 Source catalogs
tionsarepointsdrawnfromaPoissonintensitydistribution
on the celestial sphere. The intensity distribution was ob- As potential source catalogs, we consider AGNs from the
tained with a computational UHECR model. Watson et al. VCV, Swift-BAT and G10 catalogs, and galaxies from the
(2011) report strong evidence of a UHECR signal from the 2MRS catalog. This allows us to compare our analysis for
VCV AGNs. They find a low AGN fraction that is consist- the Swift-BAT and 2MRS sources with the analysis from
entwithAbreuetal.(2010).Forfiducialvaluesofthemodel Abreu et al. (2010), our analysis for the VCV sources with
parameters,theyreporta68%credibleintervalfortheAGN theanalysesfrombothAbreuetal.(2010)andWatsonetal.
fraction of FAGN =0.15+−00..1007. (2011), and our analysis of the G10 sources with Soiaporn
Soiapornetal.(2013)developedamulti-levelBayesian et al. (2013).
framework to attempt to associate the 69 UHECRs that We use the 12th edition of the VCV catalog, selecting
were recorded at the PAO in the period 2004-2009 with 17 sources with z (cid:54)0.03, as AGNs with higher redshift are
obs
nearby AGNs catalogued by Goulding et al. (2010) (here- too far away to be plausible UHECR sources, and can be
after G10). They report evidence for a small but nonzero shown to have a negligible effect on the results. We omit
fraction of the UHECRs to have originated at the AGNs sources for which absolute magnitudes are not stated. The
from G10, of the order of a few percent to 20%. total number of VCV AGNs that meet those requirements
We extend the formalism of Watson et al. (2011) with is N = 921. This is the same sample of sources that
VCV
bothagreaterdatasetandarefinedUHECRmodel.Follow- wasusedinAbrahametal.(2007),Abreuetal.(2010)and
ingAbreuetal.(2010),weextendtheanalysistotwofurther Watsonetal.(2011),andinPAO’smorerecentanalysisAab
sourcecatalogs:AGNsfromtheSwiftBurstAlertTelescope et al. (2015). While the VCV catalog is heterogenous and
(Swift-BAT) (Baumgartner et al. 2010) and galaxies from thus not ideal for statistical studies, it is close to complete
2MRS (Huchra et al. 2012). We also extend the analysis to for the low-redshift AGNs that are of relevance here.
the 17 AGNs from the G10 catalog. For the Swift-BAT catalog, we use the 58 month ver-
After discussing the UHECR and source data sets in sion, that includes a total of N = 1092 sources. In the
BAT
Section 2, we explain our UHECR model in Section 3, dis- case of the 2MRS catalog, we used the catalog version 2.4,
cuss the statistical formalism of our Bayesian model com- 2011 Dec 16. We exclude events that are within 10◦ of the
parison in Section 4, and the application of the formalism Galacticplane,toavoidbiasesduetotheincompletenessof
to mock data sets in Section 5. The results of applying the catalog in the region of the Galactic plane. This leaves
the formalism to the PAO data are discussed in Section 6. atotalofN =20,702galaxies.ThesesamplesofSwift-
2MRS
Some aspects of our computational approach are described BATand2MRSsourcesarethesameasthoseusedbyAbreu
in Appendix A, and some subtleties of our model compar- et al. (2010).
(cid:13)c 2016RAS,MNRAS000,1–14
Bayesian analysis of cosmic rays 3
0.05 0.1 0.2 0.5
Figure 1. Relative PAO exposure in Galactic coordinates. The arrival directions of the 69 UHECRs are shown as black points. The
Galacticcentre(GC)andsouthcelestialpole(SCP)areindicated.
(A)VCV (B)Swift-BAT
(C)2MRS (D)G10
0.05 0.1 0.2 0.5
Figure2.Positionaldependenceoftheexpectednumberofsourceoriginatingevents,fortheVCV,Swift-BAT,2MRS,andG10catalogs.
A fiducial value of the smearing parameter σ =3deg is assumed. The arrival directions of the 69 UHECRs are shown as black points.
Galacticcoordinatesareused,andtheGalacticcentre(GC)andsouthcelestialpole(SCP)areindicated.
TheG10catalogisawell-characterizedvolume-limited formalism(Section4),andtocreatesimulatedmockcatalogs
sample of AGNs. The 17 AGNs contained in it constitute of UHECRs to test our methods (Section 5).
all infrared-bright AGNs within 15 Mpc. This is the same
sample that was used by Soiaporn et al. (2013).
3.1 Injection
WeadoptamodelinwhichanygivenUHECRsourceemits
UHECRs with an emission spectrum given by
3 UHECR MODEL
dN /dE ∝E−γ−1, (1)
A Bayesian UHECR analysis requires a realistic model of emit emit emit
UHECR injection, propagation, and detection. This model where the logarithmic slope γ is taken to be 3.6 (Abraham
was used both to compute the likelihoods in our statistical etal.2010).Thespectrumisnormalizedinsuchawaythat
(cid:13)c 2016RAS,MNRAS000,1–14
4 A. Khanin & D. J. Mortlock
thetotalemissionrateofUHECRswithenergygreaterthan
104
E is given by
emit
dN (> E ) (cid:18)E (cid:19)−γ
emit emit =Γ emit , (2)
dt s Emin 103
where E =5.7×1019eV is the minimum UHECR emis- Ltot
min
L , Watson et al.
sion energy and Γ is the rate at which source s emits c tot
s p
UHECRs with Eemit >Emin. /M102 LLtGoZt Kfor z=0.1
L
L
adi,BH
3.2 Energy loss during propagation 101
E
thres
The energy loss processes experienced by UHECRs can
be characterized in terms of the loss length L =
loss
−E(dE/dr)−1. Given the loss length as a function of en- 1010019 1020 1021 1022
ergy, it is possible to calculate the total amount of energy E/eV
that a UHECR loses as it travels to the Earth from a given
distance by solving the differential equation Figure3.Losslengthsfromthethreeenergylossprocesses,com-
paredtotheconstantconstantlosslengthusedbyWatsonetal.
dE E
=− . (3) (2011),asdescribedinSection3.2.
dr L (E)
loss
For pure proton composition, L obeys the expression
loss 3.3 Effective smearing
1
L−los1s = c[βadi(E,z)+βGZK(E,z)+βBH(E,z)], (4) WecombinedthemagneticdeflectionthataUHECRexperi-
encesduringpropagationandtheuncertaintyinitsdetected
wherecisthespeedoflightandβGZK(E,z),βBH(E,z)and arrivaldirectionintoasinglekernel,whichwaschosentobe
βadi(E,z)aretermscorrespondingtothethreemainenergy a von Mises-Fisher (vMF) distribution, defined as
lossprocessesexperiencedbyUHECRsofpureprotoncom-
κ
position (e.g. Stanev 2009): Pr(rˆ|rˆsrc,κ)= 4πsinh(κ)exp(κrˆ·rˆsrc), (5)
(i) the GZK scattering off the CMB photons at energies where rˆis the measured arrival direction of the ray, rˆ is
src
above E (cid:38)5×1019eV; the source direction and κ is the concentration parameter.
(ii) Bethe-Heitler(BH)e+e−pairproduction(alsoascat- The vMF distribution resembles a Gaussian on the sphere,
tering process off the CMB radiation), which dominates at withκbeinginverselyrelatedtothewidthoftheGaussian:
lower energies (Hillas 1968); forlargevaluesofκthedistributionispeakedoveranangu-
√
(iii) theadiabaticenergylossduetotheexpansionofthe larscaleof∼1/ κ;ifκtendsto0thedistributionbecomes
Universe. uniform on the sphere.
The magnitude of the deflection that the highest en-
A detailed discussion of these terms, including expressions
ergyUHECRsexperienceisuncertain,withtheestimatesof
andparametrizations,canbefoundinDeDomenico&Inso-
typical deflection angles ranging from ∼2 to ∼10 deg (e.g.
lia(2013).Fortheenergiesthatarerelevantinthisinvestig-
MedinaTancoetal.1998;Sigletal.2004;Dolagetal.2005).
ation, the dominant term is β (E,z). The Bethe-Heitler
GZK Weassumeafiducialsmearingangleofσ(cid:39)3deg(κ=360),
and adiabatic processes dominate the energy loss at lower
but also conduct investigations for smearing angles of σ (cid:39)
energies, but play only a minor role at the higher energies
6 and 10 deg (κ=90 and 30).
in question.
The loss lengths are shown as a function of energy in
Figure 3. The contributions to the loss length from the BH 3.4 Observed UHECR flux
and adiabatic losses are combined into a single function
The number of UHECRs from source s above a threshold
L that is contrasted with the loss length due to the
adi,BH
energyE observedonEarthperunitareaperunittime,
GZKeffect,L .Thetwoarecombinedintothetotalloss thres
GZK dN (E (cid:62)E )/dtdA,isaquantitythatisimportantin
length L . The figure shows L plots for z values of 0.0 s obs thres
tot tot
ourstatisticalanalysis.Thisrateisproportionaltotherate
and0.1,whichcorrespondtodistancesof0and∼400Mpc,
of UHECRs emitted by the source, Γ , but it also depends
thus covering the GZK horizon. L appears very rapidly s
GZK
after an energy of ∼ 4×1019eV and begins to dominate onthedistance-dependenceoftheUHECRenergyloss,and
on the UHECR injection spectrum. We use the UHECR
the energy loss. As we are interested only in UHECRs with
energies E >E =5.7×1019eV, the GZK scattering propagation model described in Section 3.2 to determine
obs thres
the injection energy corresponding to the threshold energy
is the most relevant loss process in this investigation.
E and to the source distance D . Combining this value
The energy dependence of L (E) is one of the main thres s
loss
withEquation2andwiththesourcedistanceD ,weobtain
improvementsofthispropagationmodeloverthemodelused s
in Watson et al. (2011), where L was taken to be a con-
loss
stant. The constant value of Lloss used by Watson et al. dNs(Eobs (cid:62) Ethres) = Γs (cid:20)Eemit(Ethres)(cid:21)−γ. (6)
(2011) is also displayed in Figure 3 for comparison. dtdA 4πD2 E
s min
(cid:13)c 2016RAS,MNRAS000,1–14
Bayesian analysis of cosmic rays 5
This expression assumes that the observed energy E is 4.2 The likelihood
obs
equivalenttothearrivalenergyoftheUHECR,E .Thus,
arr To compute the likelihood, we use a ‘counts in cells’ ap-
for the purposes of the calculation, the 12% energy uncer-
proach, in which the sky is divided into 1800 × 3600 =
taintyofthePAOmeasurementsisneglected.Thevariation
6,480,000 pixels, that are distributed uniformly in right as-
in source rates Γ among the sources that we are consider-
s censionanddeclination.Thus,thedatasetd canberewrit-
ing is not negligible. We use the source rate of Centaurus
ten as a set of counts in each pixel {N }.
A as the reference value Γ. The source rate of a source s is c,p
ThelikelihoodPr(d|Γ,R)isthengivenbyaproductof
obtained by weighing the flux F of that source in a par-
s the individual Poisson likelihoods in each pixel, and can be
ticular band against the flux F of Centaurus A in that
Cen written as
same band. The wave band of the flux thereby is different
depending on the source catalog. For VCV, the flux of the Pr(d|Γ,R)
sourceintheV-bandisused,forSwift-BATtheX-rayflux,
for 2MRS the IR flux, for G10 the K-band flux. The fluxes = (cid:89)Np (Nsrc,p+Nbkg,p)Nc,pexp[−(Nsrc,p+Nbkg,p)], (10)
are thus used as weights, so that sources with higher flux N !
c,p
contribute more UHECRs. This approach is very similar to p=1
theapproachusedinAbreuetal.(2010),wherefluxeswere
where N and N are the expected counts in pixel p
src,p bkg,p
used to weigh the sources from the Swift-BAT and 2MRS
due to sources and background, respectively. The expected
catalogs in the same way. Incorporating the fluxes into the
number of counts in pixel p that are contributed by the
formalism, we obtain the expression
background is
dNs(Eobs (cid:62) Ethres) = Γ Fs (cid:20)Eemit(Ethres)(cid:21),−γ (7) N =R(cid:90) d(cid:15) dΩ , (11)
dtdA 4πDC2enFCen Emin bkg,p p dΩ obs
where DCen is the distance to Centaurus A. wheretheintegralisoverthepixelp,andd(cid:15)/dΩistherelat-
ive exposure (Section 2.1). The expected number of source
originating events in pixel p is
4 STATISTICAL FORMALISM N =(cid:88)Ns dNs(Eobs (cid:62)Ethres)(cid:90) d(cid:15)Pr(˜r |˜r )dΩ , (12)
src,p dtdA dΩ obs s obs
Given a sample of UHECRs arrival directions, we would s=1 p
like to determine the fraction of these rays that have come where the sum is over the sources, Pr(˜r |˜r ) is the vMF
obs s
from a set of sources under consideration. To do so, we use distribution (Equation 5), and dN (E (cid:62) E )/dtdA is
s obs thres
a two-component parametric model characterized by two theobservedUHECRfluxdiscussedinSection3.4.Inserting
rates: The source rate Γ and the isotropic background rate Equations 11 and 12 into Equation 10, we arrive at the full
R. As elaborated in Section 3.4, we use the source rate of likelihood.
Centaurus A as the reference value of Γ. We obtain a joint ThepositionaldependenceofN followstherelative
bkg,p
posterior distribution for the two rates: exposure of PAO, as shown in Figure 1. The positional de-
pendence of N depends both on the PAO exposure and
Pr(Γ,R)Pr(d|Γ,R) src,p
Pr(Γ,R|d)= (cid:82)∞(cid:82)∞Pr(Γ,R)Pr(d|Γ,R)dΓdR, (8) onthedistributionofsourcesinthegivencatalog.Figure2
0 0 showsthedependenceforthefourcatalogsthatareusedin
where Pr(Γ,R) is the prior distribution for Γ and R, and thisstudy.Thedependenceisdominatedbythedistribution
Pr(d|Γ,R)isthelikelihood(i.e.theprobabilityofobtaining of local AGNs, by far the strongest source being Centaurus
the data set d given values of Γ and R). A(l=309.5◦,b=19.4◦),whichpreviouslystudies(e.g.Ab-
rahametal.2007)havesuggestedasthedominantUHECR
source.
The expression for the likelihood can be rearranged to
4.1 Prior
reduce the total number of computations, as described in
We adopt a uniform prior over Γ and R, with Γ(cid:62)0, R(cid:62)0. Appendix A.
This plausibly encodes our ignorance of the two paramet-
ers,and,unlikemaximumentropypriors,includesapossible
value of 0 for both parameters. The maximum values of Γ 4.3 The source fraction
and R are denoted as Γ and R . We have conducted
max max The source fraction1 is defined as the fraction of the
our analysis for flat priors of varying width, using a vari-
UHECRsthatareexpectedtohaveoriginatedatthesources
ablewidthparameters.Theexpressionforthepriorcanbe
inwhichevercatalogisunderconsiderationandisgivenby
written as
Pr(Γ,R|d,M )= 1 . (9) F (Γ,R)= (cid:80)Np=p1Nsrc,p . (13)
2 s2ΓmaxRmax src (cid:80)Np=p1Nsrc,p+Nbkg.p
Γ and R have been chosen in such a way that when
max max
s=1,thepriorcoversthe99.7%credibleregionimpliedby
the likelihood and an infinitely broad uniform prior. This 1 The source fraction Fsrc is equivalent to the AGN fraction
gives a data driven scaling for the rates. The priors and FAGN used in Watson et al. (2011) but now generalized to al-
their dependence on s are illustrated in Appendix B. lowfornon-AGNprogenitors.
(cid:13)c 2016RAS,MNRAS000,1–14
6 A. Khanin & D. J. Mortlock
The posterior for F can be calculated from the posterior set of 69 events should be sufficient to distinguish between
src
over the rates as the two models. Figure 4 also shows the Bayes factors as
functions of s for the two cases. The Bayes factors B
Pr(F |d) 21
src that are displayed are the inverses of the SDDR given in
Γ(cid:90)maxR(cid:90)max Equation17,andfavourthemorecomplexmodelforBayes
= Pr(Γ,R|d)δ [F −F (Γ,R)]dΓdR. (14) factors >1.
D src src
To assess the results of the Bayes factor simulations,
0 0
we can derive a rough range of plausible values of s from
Pr(Fsrc|d)isinsensitivetoRmaxandΓmaxprovidedtheyare physical models, and then look at the behaviour of the
sufficiently large. Bayes factors at those physically plausible values. Plaus-
ible models of UHECR injection predict that the UHECR
luminosity of a source like Centaurus A is of the order of
4.4 Model comparison
2.9×1039ergs−1 (cid:39)1.81×1051eV(Fraijaetal.2012).Ifthis
WewouldliketocomparemodelM wherealltheUHECRs istakenasthetypicalUHECRluminosityofasource,then
1
aredrawnfromauniformdistributionwithmodelM where foraUHECRenergyrangeof(5.7−100)×1019eV,therange
2
the UHECRs are derived from a combination of a back- of source rates can be calculated by dividing the UHECR
ground and a source originating component. To do this, we luminosity by the limiting values of this range. The result
conductaBayesianmodelcomparison.Foradatasetd,and of this calculation is a range of source rates Γ of roughly
twomodelsM andM ,theratioofthemarginallikelihoods (2−33)×1030s−1. The values of s corresponding to this
1 2
for the two models, termed the Bayes factor, is rangehavebeenmarkedonFigure4.(Thevaluesareslightly
different for each of the simulations. For the sake of clarity,
Pr(d|M )
B = 1 . (15) onlythevaluesfortheuniformsimulationaredisplayed,the
12 Pr(d|M )
2 others being broadly similar.) For the sourced case, model
In the specific case that is considered here, the models M isstronglyfavouredforallphysicallyplausiblevaluesof
2
arenested:WhenΓ=0,modelM reducestomodelM .A s,whilefortheuniformcase,thesimpleuniformmodelM
2 1 1
general expression of the Bayes factor in this situation is is favoured for the physically plausible values.
(cid:82)
Pr(R|M )Pr(d|R,M )dR
B12 = (cid:82) Pr(Γ,R|M )1Pr(d|Γ,R,M1 )dΓdR. (16)
2 2
Itcanbeshown(Dickey1971)thatinthecaseofsuchnested
6 RESULTS
models, the expression reduces to
Theresultsoftheapplicationofthestatisticalmethodsde-
Pr(Γ=0|d,M )
B = 2 . (17) scribed in Section 4 to the data described in Section 2 are
12 Pr(Γ=0|M )
2 showninFigures5and6.Figure5contraststheresultsfrom
This expression is known as the Savage-Dickey Density Ra- our analysis with the equivalent results from Watson et al.
tio,orSDDR.Qualitatively,thisexpressionmeansthatthe (2011), and with the results for an intermediate case. The
nested uniform model is preferred if, within the context of use of a more refined propagation model leads to a higher
themorecomplicatedmodel,thedataresultinanincreased posterior probability for lower source rates. The reason for
probability that Γ=0. thatisthatinWatson’spropagationmodel,theenergyloss
length is constant and very small (Figure 3). UHECRs ex-
perience more drastic energy loss than in the more realistic
model,whichleadstomoredistantAGNsbeingexcludedas
5 SIMULATIONS
plausible source candidates. As fewer sources are included,
In order to investigate the constraining power of a data set ahighersourcerateisrequiredtogeneratethesamesample
of 69 events, we apply the method to simulated data sets. of UHECRs.
We use two extreme cases: The inclusion of 69 events reduces the extent to which
thenon-uniformmodelisfavoured.Thisisevidentfromthe
(i) Uniform arrival directions. These rays were drawn
posterior of the source fraction, and also from the behava-
from a probability distribution that followed the PAO ex-
viour of B . This result agrees with the results of Abreu
posure. 21
et al. (2010), which reported that the full 69 events yield
(ii) UHECRs originating at sources from a catalog. We
lower evidence of anisotropy than the earlier study Abra-
conducted simulations for all four of the catalogs. In each
ham et al. (2007), which analysed 27 events.
catalog, the sources were weighted by their fluxes and the
Figures 6 shows results for all four of the source cata-
PAOexposure.Randomsourceswerethenselected,andthe
logs,andforallvaluesofthesmearingparameter.Displayed
propagation model of Section 3 was used to propagate rays
are the posteriors for the source fraction, as well as plots of
from the sources to the Earth.
B againsts.Theconstraintsonthesourcefractionforall
21
The posteriors for the source and background rates, as casesareshowninTable1.Thefiguresandtableshowthat
well as the posteriors for the source fraction, are summar- for greater smearing, the range of plausible values of F is
src
ized in Figure 4. The posteriors for the uniform and source increased,andthemostprobablevalueofthesourcefraction
centred cases are completely disjoint, which demonstrates ishigherthanforthefiducialmodelofσ=3deg.Thereason
thatinextremesenarioswhereallUHECRsoriginateeither is that for greater magnetic deflection, the UHECR intens-
fromauniformbackgroundorfromasourcecatalog,adata ity distribution becomes more uniform, so that the uniform
(cid:13)c 2016RAS,MNRAS000,1–14
Bayesian analysis of cosmic rays 7
Table 1. Maximum a posteriori estimates and 68% credible in-
(A)
tervalsforFsrc.
Uniform simulation
VCV simulation
30)060 Swift-BAT simulation Catalog σ=3deg σ=6deg σ=10deg
1× 2MRS simulation
1− G10 simulation VCV 0.09+−00..0054 0.14+−00..0076 0.22+−00..0098
s
1− Swift-BAT 0.25+0.09 0.37+0.11 0.46+0.13
rc40 −0.08 −0.10 −0.12
(s 2MRS 0.24+0.12 0.33+0.14 0.40+0.15
Γ −0.10 −0.14 −0.15
ate G10 0.08+−00..0043 0.14+−00..0065 0.22+−00..0077
r
ce 20
r
u
o
s
0
0 10 20 30 40
background rate R(sr−1m−2s−1 10−17)
×
andmixedmodelsbecomemoredifficulttodistinguish,and
(B)
20 a greater range of F values become viable.
src
Uniform simulation
The plots of B demonstrate that for all physically
VCV simulation 21
plausible prior ranges of the model parameters, the fully
Swift-BAT simulation
15 2MRS simulation isotropic model is disfavoured. The form of the dependence
G10 simulation of B21 on s is elaborated upon in Appendix B.
TheseresultsfortheVCV,Swift-BAT,and2MRScata-
)
Fdsrc|10 lwohgsocuasnedbeaccoomrrpealraetdiown-ibthastehdearnesaulyltssisofoAnbtrheeuVetCaVl.c(a2t0a1l0o)g,
(
Pr that mirrored the analysis in Abraham et al. (2007). Ab-
reu et al. (2010) reported a correlation of (38+7)% between
−6
UHECRsandsourcesfromtheVCVcatalog,whichwascon-
5
siderablylowerthanthanthe(69+11)%correlationthatwas
−13
reported in Abraham et al. (2007). This reduction in the
correlation is consistent with our findings that the source
00.0 0.2 0.4 0.6 0.8 1.0 fraction is reduced as we increase the data set from 27 to
Fsrc 69 events. In addition to these correlation based methods,
Abreuetal.(2010)conductedalikelihoodbasedstudysim-
(C) ilartotheanalysispresentedhere,wherethelikelihoodwas
1045 Uniform simulation takenasaprobabilitymapofarrivaldirectionsofUHECRs,
1041 VCV simulation parametrized by a magnetic smoothing angle σ and a frac-
1037 Swift-BAT simulation tion of isotropic rays fiso, which is equivalent to 1−Fsrc.
1033 2MRS simulation Theselikelihood-basedstudieswereconductedfortheSwift-
G10 simulation
1029 BATand2MRScatalogs.Forthe2MRScase,themaximum
1025 likelihoodvaluesoffiso andσ arereportedas0.56and7.8◦,
B211021 respectively. The σ value lies between our chosen smearing
1017 angles 6deg and 10deg. The value for fiso corresponds to
a value of F of 0.44, which is consistent with our F
1013 src src
credible intervals for these chosen smearing angles. For the
109
case of Swift-BAT, the maximum likelihood value of f is
105 iso
givenas0.64,whichcorrespondstoasourcefractionof0.36.
101 M2 favoured
M1 favoured The maximum likelihood estimate of the smearing angle is
10-3
reported as 1.5◦, which is lower than our minimum chosen
10-7
10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 value of 3deg. Despite the difference between the angles, a
s
F value of 0.36 can still be considered broadly consistent
src
Figure 4. Results from simulations: Uniform UHECRs, and with the 68% credible interval for 3deg, 0.25+0.09.
−0.08
UHECRsoriginatingatsourcesfromthe4catalogs.Inallcases, Our results for the G10 catalog can be compared with
69eventsareused.(A)PosteriorsforΓandR.Thecontoursare theworkofSoiapornetal.(2013).Thatanalysisinvolvedthe
the 68.3%, 95.4% and 99.7% highest posterior density credible
full data set of 69 events, and found evidence for small but
regions.(B)Posteriorsforthesourcefraction.(C)PlotofBayes
nonzerovaluesofF ,oftheorderofafewpercentto20%,
factorsB21 asafunctionofthehyperparameters.In(C),the×- src
ruling out values of F > 0.3. This is broadly consistent
markandtheverticallinesignifytheminimumandthemaximum src
withourresults,whichsuggestthatvaluesofF <∼0.3are
values of the physically plausible range of s. The minimum and src
maximum values that are displayed correspond to the uniform the most probable for all values of the smearing parameter.
simulation.
(cid:13)c 2016RAS,MNRAS000,1–14
8 A. Khanin & D. J. Mortlock
7 CONCLUSIONS
(A)
We have performed a Bayesian analysis of the 69 UHECRs
69 events, variable L
loss detected by the PAO with energies E > 5.7×1019eV
27 events, variable L obs
) loss to determine the fraction of these UHECRs that originated
3010× 27 events, constant Lloss fromcatalogsofplausibleUHECRsources.Thesourcescon-
sidered were AGNs from the VCV, Swift-BAT, and G10
1−
s catalogs, and galaxies from the 2MRS catalog.
1−
rc20 For the fiducial magnetic smearing parameter of σ= 3
s
( deg,wereport68%credibleintervalsforthesourcefraction
Γ
e of0.09+0.05,0.25+0.09,0.08+0.04 and0.24+0.12 fortheVCV,
t −0.04 −0.08 −0.03 −0.10
ra Swift-BAT, G10 and 2MRS catalogs, respectively. For all
ce physicallyplausiblevaluesofthemodelparameters,thefully
r
u uniform model is disfavoured. The results of our study are
o
s inbroadagreementwithpreviousworkonthissubject,such
as Watson et al. (2011), Abreu et al. (2010) and Soiaporn
00 10 20 et al. (2013). The credible intervals for the VCV catalog
background rate R(sr−1m−2s−1 10−17) arelowerthantheanalogouscredibleintervalsfromWatson
×
etal.(2011),whichusedasimilarmethodtoanalyse27PAO
(B) events. This is consistent with earlier studies: Abreu et al.
8
F =0.09+0.05 69 events, variable Lloss (2010), which analysed 69 events, reported a lower signal
7 src −0.04 27 events, variable Lloss of anisotropy than the earlier study Abraham et al. (2007),
27 events, constant Lloss which used 27 events.
6 We will extend this Bayesian framework to include the
arrival energies of the UHECRs as well as the arrival direc-
Fd)src|45 Fsrc=0.15−+00..0170 tionsI.tisexpectedthatfutureexperimentswillproducedata
( sets that will be sufficiently large for our Bayesian method
r
P
3 (and other statistical approaches; see e.g. Rouill´e d’Orfeuil
etal.2014)todetecteventheweakclusteringexpectedifthe
2 UHECRShavecomefromnearbysources.PAOiscontinuing
to take data and is expected to produce a sample of ∼250
1
UHECRsoveritsfirstdecadeofoperations.Lookingfurther
ahead, the planned Japanese Experiment Module Extreme
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Universe Space Observatory (JEM-EUSO, Adams Jr. et al.
F
src
2013) on the International Space Station (ISS) is scheduled
forlaunchin2017andisexpectedtodetect∼200UHECRs
(C)
109 annually over its five year lifetime.
108
107 References
106
Aab A., et al., 2015, The Astrophysical Journal, 804, 15
105
AbbasiR.,etal.,2010,TheAstrophysicalJournalLetters,
104
713, L64
103
B21102 Abbasi R. U., et al., 2008, Physical Review Letters, 100,
101101
101
100 MM21 ffaavvoouurreedd Abbasi R. U., et al., 2008, Astroparticle Physics, 30, 175
Abraham J., Abreu P., Aglietta M., Ahn E. J., Allard D.,
10-1
AllenJ.,Alvarez-Mun˜izJ.,AmbrosioM.,AnchordoquiL.,
10-2
10-3 69 events, variable Lloss Andringa S., et al. 2010, Physics Letters B, 685, 239
27 events, variable Lloss Abraham J., et al., 2007, Science, 318, 938
10-4 27 events, constant Lloss Abraham J., et al., 2008, Astroparticle Physics, 29, 188
10-5
10-4 10-3 10-2 10-1 100 101 102 103 104 Abraham J., et al., 2008, Physical Review Letters, 101,
s
061101
Figure 5.Resultsforσ=3deg,andthesourcesfromtheVCV Abraham J., et al., 2009, arxiv:0906.2347
catalog.Resultsfor27and69events,andforconstantandvari- Abreu P., et al., 2010, Astroparticle Physics, 34, 314
ablelosslengthsaredisplayed.(A)Posteriorsforthesourceand Abu-Zayyad T., et al., 2012, The Astrophysical Journal,
backgroundrates.Thecontoursarethe68.3%,95.4%and99.7%
757, 26
highest posterior density credible regions. (B) Posterior for the
Adams Jr. J. H., et al., 2013, arXiv:1307.7071
source fraction. (C) Plot of Bayes factors B21 as a function of
Baumgartner W. H., Tueller J., Markwardt C., Skinner
thehyperparameters.In(C),physicallyplausiblerangesofsare
G.,2010,inAAS/HighEnergyAstrophysicsDivision#11
shown for the cases of 27 events (blue) and 69 events (black),
with a variable loss length. The ×-marks and the vertical lines Vol.42ofBulletinoftheAmericanAstronomicalSociety,
signify the minimum and the maximum values of the physically The Swift-BAT 58 Month Survey. p. 675
plausiblerangesofs.
(cid:13)c 2016RAS,MNRAS000,1–14
Bayesian analysis of cosmic rays 9
De Domenico M., Insolia A., 2013, Journal of Physics G
Nuclear Physics, 40, 015201
DolagK.,GrassoD.,SpringelV.,TkachevI.,2005,J.Cos-
mology & Astro-Part. Phys., 1, 9
Fraija N., et al., 2012, The Astrophysical Journal, 753, 40
George M. R., Fabian A. C., Baumgartner W. H.,
Mushotzky R. F., Tueller J., 2008, MNRAS, 388, L59
GouldingA.D.,AlexanderD.M.,LehmerB.D.,Mullaney
J. R., 2010, MNRAS, 406, 597
GregoryP.C.,2010,BayesianLogicalDataAnalysisforthe
Physical Sciences. Cambridge University Press, pp 376–
388
Greisen K., 1966, Physical Review Letters, 16, 748
Hillas A. M., 1968, Canadian Journal of Physics, 46, 623
HuchraJ.P.,etal.,2012,TheAstrophysicalJournal,199,
26
Ivanov A. A., 2009, Nuclear Physics B Proceedings Sup-
plements, 190, 204
Kalmykov N. N., Khrenov B. A., Kulikov G. V., Zotov
M. Y., 2013, Journal of Physics Conference Series, 409,
012100
Kotera K., Olinto A. V., 2011, ARA& A, 49, 119
Letessier-SelvonA.,etal.,2014,BrazilianJournalofPhys-
ics, 44, 560
Letessier-Selvon A., Stanev T., 2011, Reviews of Modern
Physics, 83, 907
MedinaTancoG.A.,deGouveiaDalPinoE.M.,Horvath
J. E., 1998, The Astrophysical Journal, 492, 200
Pe’ErA.,MuraseK.,M´esza´rosP.,2009,Phys.Rev.D,80,
123018
Rouill´e d’Orfeuil B., Allard D., Lachaud C., Parizot E.,
Blaksley C., Nagataki S., 2014, arXiv:1401.1119
SiglG.,MiniatiF.,EnßlinT.A.,2004,PhysicalReviewD,
70, 043007
SoiapornK.,ChernoffD.,LoredoT.,RuppertD.,Wasser-
man I., 2013, The Annals of Applied Statistics, 7, 1249
Stanev T., 2009, New Journal of Physics, 11, 13
Tinyakov P. G., Tkachev I. I., 2001, Soviet Journal of Ex-
perimental and Theoretical Physics Letters, 74, 445
V´eron-CettyM.-P.,V´eronP.,2006,AstronomyandAstro-
physics, 455, 773
Watson L. J., Mortlock D. J., Jaffe A. H., 2011, MNRAS,
418, 206
Zatsepin G., Kuzmin V., 1966, JETP Lett., 4
(cid:13)c 2016RAS,MNRAS000,1–14
10 A. Khanin & D. J. Mortlock
(A) VCV (A) VCV
87 Fsrc=0.09−+00..0045 σσ == 36 ddeegg 11001102 σσ == 36 ddeegg
6 Fsrc=0.14−+00..0067 σ = 10 deg 108 σ = 10 deg
Fd)src|45 Fsrc=0.22−+00..0089 B21 110046
Pr(3 102
2 100 MM21 ffaavvoouurreedd
1 10-2
0 10-4
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 10-4 10-3 10-2 10-1 100 101 102 103
F s
src
(B) Swift-BAT (B) Swift-BAT
5
F =0.25+0.09 σ = 3 deg 1022 σ = 3 deg
src −0.08 σ = 6 deg 1020 σ = 6 deg
4 Fsrc=0.37+00..1101 σ = 10 deg 11001168 σ = 10 deg
Fd)src|3 Fsr−c=0.46−+00..1123 B21111000111024
Pr(2 110068
104
1 102
100 MM21 ffaavvoouurreedd
0 10-2
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 10-4 10-3 10-2 10-1 100 101 102 103 104 105
F s
src
(C) 2MRS (C) 2MRS
323...055 Fsrc=F0s.r2c4=−+000...113023+00..1144 σσσ === 361 0dd eedggeg 1111000011124680 σσσ === 361 0dd eedggeg
Fd)src|2.0 Fs−rc=0.40−+00..1155 B2111100011802
Pr(1.5 110046
1.0 102
0.5 1100-02 MM21 ffaavvoouurreedd
0.0 10-4
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 10-4 10-3 10-2 10-1 100 101 102 103
Fsrc s
(D) G10 (D) G10
10 F =0.08+0.04 σ = 3 deg 1014 σ = 3 deg
src −0.03 σ = 6 deg 1012 σ = 6 deg
8
F =0.14+0.06 σ = 10 deg 1010 σ = 10 deg
src 0.05
Fd)src|6 Fsrc=−0.22+00..0077 B21 110068
Pr( 4 − 104
102
2 100 MM21 ffaavvoouurreedd
10-2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 10-140-4 10-3 10-2 10-1 100 101 102 103
Fsrc s
Figure 6. Posteriors of the source fraction, and plots of B21 against the hyperparameter s, for the three smearing angles σ = 3deg,
6degand10deg,andforthethreesourcecatalogs(A)VCV,(B)Swift-BAT,and(C)2MRS.TheplotsofB21showphysicallyplausible
rangesofs:The×-marksandtheverticallinessignifytheminimumandthemaximumvaluesoftheseranges.
(cid:13)c 2016RAS,MNRAS000,1–14