ApJ, 679, 587, 2008

20 Pages·2008·12.53 MB·English

TheAstrophysicalJournal,679:587Y606,2008May20 A
#2008.TheAmericanAstronomicalSociety.Allrightsreserved.PrintedinU.S.A.
LINE SEARCHES IN SWIFT X-RAY SPECTRA
C. P. Hurkett,1 S. Vaughan,1 J. P. Osborne,1 P. T. O’Brien,1 K. L. Page,1 A. Beardmore,1 O. Godet,1
D. N. Burrows,2 M. Capalbi,3 P. Evans,1 N. Gehrels,4 M. R. Goad,1 J. E. Hill,4,5
J. Kennea,2 T. Mineo,6 M. Perri,3 and R. Starling1
Received2007April13;accepted2008January29
ABSTRACT
Priortothelaunchof theSwift missionseveralX-raylinedetectionswerereportedingamma-rayburstafterglow
spectra.Todate,thesepreYSwifteraresultshavenotbeenconclusivelyconfirmed.Themostcontentiousissueinthis
areaisthechoiceof statisticalmethodusedtoevaluatethesignificanceof thesefeatures.Inthispaperwecompare
threedifferentmethodsalreadyextantintheliteratureforassessingthesignificanceofpossiblelinefeaturesanddis-
cusstheirrelativeadvantagesanddisadvantages.Themethodsaredemonstratedbyapplicationtoobservationsof
40burstsfromthearchiveofSwiftXRTatearlytimes(<afewkilosecondsposttriggerintherestframeoftheburst).
Basedonthisthoroughanalysiswefoundnostrongevidenceforemissionlines.Foreachof thethreemethodswe
havedetermineddetectionlimitsforemission-linestrengthsinburstswithspectralparameterstypicaloftheSwift-era
sample.Wealsodiscusstheeffectsofthecurrentcalibrationstatusonemission-linedetection.
Subject headingsg: gamma rays: bursts — gamma rays: observations — line: identification —
methods: data analysis — methods: statistical
Online material: color figures
1. INTRODUCTION The models for the production of such emission features
are divided into transmission and reflection models, although
It is widely accepted that the spectra of the X-ray afterglow
thelargeequivalentwidths((cid:1)afewkeV)inferredfromtheob-
ofgamma-raybursts(GRBs)aredominatedbynonthermalemis-
servedX-rayfeaturesfavormodelsinwhichthelineisproduced
sion, the leading candidate for which is synchrotron emission
by reflection (Rees & Me´sza´ros 2000; Ballantyne & Ramirez-
(Piran 2005,andreferencestherein), although alternateemis-
Ruiz2001;Vietrietal.2001).Proposedmodelshavetoovercome
sionprocesseshavealsobeensuggested,suchasself-Compton
two constraints; the size problem and the kinematic problem.
(Waxman1997;Ghisellini&Celotti1999)orinverseCompton
Observing a line at a time t after the burst implies that the
scatteringof externallight(Brainerdetal.1994;Shemi1994; obs
emitting material must be within a distance of (cid:1)ct /(1þz)
Shaviv&Dar1995;Lazzatietal.2004). obs
fromthecentralengine,thusimplyingthattheregionmustbe
Up to thepresenttime theX-rayspectra of Swift afterglows
compactif aline,orlines,areobservedatearlytimes.Inaddi-
aregenerallywelldescribedbyanabsorbedpowerlaw(forcoun-
tion, the emitting region must contain (cid:1)0.1 M of Fe (in the
terexamplesseeButler2007),typicallyabsorbedbymaterialwith (cid:2)
caseof FeK(cid:1)features)whilestillbeingopticallythintoelec-
acolumndensityinexcessof thewell-measuredGalacticvalues
tronscattering,inorderthatComptonizationdoesnotbroaden
(Campanaetal.2006c).Table2ofCampanaetal.(2006c)shows
thelinebeyondtheobservedwidths(Vietrietal.2001).If the
that,of17burstsanalyzed,14haveobservedN valuesgreater
H linewidthisinterpretedasbeingduetothevelocityof thesuper-
thanthemeasuredGalacticcolumndensity,whiletheremain-
novaremnant,theobservedlimitonthiswidthimpliesanagelimit
ing three have observed N values that are consistent, within
H ontheremnantof (cid:1)10Y20days.However,atthistimeConuclei
limits,withthemeasuredvalues.
outnumberbothNiandFenuclei;thustheemissionlinewouldbe
Inthepastithasbeenproposedthatthereareotherspectralfea-
duetoCoatanenergyof 7:5/(1þz)keV,whichisthekinematic
tures,withvaryinglevelsof significance,inadditiontothebasic
problem.
absorbedpowerlawspectrum(Piroetal.1999,2000;Yoshidaetal.
Variousgeometrieshavebeensuggestedforthereflectionmod-
1999;Amatietal.2000;Antonellietal.2000;Reevesetal.2002;
els,whichrelyoneitheraprecursororsimultaneoussupernova
Watsonetal.2002,2003;Fronteraetal.2004).Mostareattributed
(SN)event.If aSNoccursseveraltensof daysbeforetheGRB
toFeK(cid:1)emissionlinesortheradiativerecombinationcontinuum
thissolvesboththesizeandkinematicproblems.Inthesecases
of thesameelement.SomehavebeenattributedtotheK(cid:1)linesof
the radiation from the GRB jets can either illuminate the inner
Ni,Co,orlighterelementssuchasSi,S,Ar,andCa.Intwocases
faceof theSNshellremnantortheinnerfacesof widefunnels
therehavebeenareportsofatransientabsorptionfeaturealsocor-
thattheyexcavatethroughyoungplerionicremnants.However,
respondingtoFeK(cid:1)(Amatietal.2000;Fronteraetal.2004).
thesemodelshavebeenquestionedfollowingthesimultaneous
GRB-SN association indicated byGRB 980425 (Galama etal.
1998)andconfirmedbyGRB030329(Hjorthetal.2003;Stanek
1 XROA Group, Department of Physics and Astronomy, University of
etal.2003)andGRB060218(Campanaetal.2006b;Pianetal.
Leicester,LeicesterLE17RH,UK;[email protected].
2 DepartmentofAstronomyandAstrophysics,PennsylvaniaStateUniversity. 2006).Inthiscasethemostlikelyscenarioforemission-linepro-
3 ASIScienceDataCenter,ASDCcareofESA-ESRIN,viaG.Galilei00044 ductionoccursif theprogenitorejectsalargeamountof matter
Frascati,Italy. atsubrelativisticspeedsalongitsequator.Thehaloof material
4 NASAGoddardSpaceFlightCenter,Greenbelt,MD20771.
surroundingmassivestars,ejectedbytheirstrongstellarwinds
5 Universities Space Research Association, 10211 Wincopin Circle, Suite
towardtheendoftheirmainsequencelifetime,scattersafraction
500Columbia,MD21044.
6 INAFIASF-Pa,viaU.LaMalfa153,90146Palermo,Italy. of the photons from the prompt and afterglow phase back into
587
588 HURKETT ETAL. Vol. 679
theequatorialmaterial,whichthenproducesX-raylineemission of 153bursts,40of whichcontainedsufficientWTmodedata
(Vietrietal.2001). forouranalysismethods.WTmodedatawaschosenprimarily
Verifyingthepresenceof suchspectralfeaturesisof critical because the time interval covered by these observations, typi-
importanceastheywillallowustoprobethecircumburstenvi- callyT þ0stoT þ500s(althoughforbrightburststhismay
ronmentof theGRB,aswellasgaininganindirectindicationof extenduptoT+fewks)intherestframeoftheburst,israrely
thepossiblestructureandbehaviorof thecentralengine. explored.Priorobservationsofthe1999Y2003burststypically
Thestatisticalsignificanceof the1999Y2003reportedfeatures startat20+hrafterthetriggerintheobserver’sreferenceframe,
islow(usually2Y3(cid:2)),andonlytwodetectionshaveasignificance althoughAntonellietal.(2000)reportonanemission-linedetec-
>4(cid:2)(GRB991216:4.7(cid:2),Piroetal.2000;andGRB030227:4.4(cid:2), tionatTþ(cid:1)12hr;Amatietal.(2000)andFronteraetal.(2004)
Watsonetal.2003).Althoughlaterdetectionsweremadewith reportabsorptionlinefeaturesinveryearlytimedata(Tþ<20s
much more sensitive instruments than the early ones, all detec- and Tþ<300 s, respectively) from the Wide Field Camera
tionshaveremainedatthislowsignificancelevel;asaresultthey (WFC)of BeppoSAX.InadditionWTmodedataisonlytaken
remainthesubjectof muchdebate.Arguablythemostcontrover- whiletheGRBafterglowisbright.Allofthemethodsdiscussed
sialissueinthediscussionoflinedetectionsisthechoiceofsta- caneasilybeextendedtophotoncounting(PC)modedata.We
tisticalmethodemployedtogaugetheirsignificance.Atleastfour acknowledgethatthecurrenttheoreticalmodelsforlineemission
methodshavebeenproposedandusedintheGRBliterature. indicatethatlinescouldoccurattimesnotcoveredbyWTmode
data,however,thesamemodelsdonotruleoutthistimeperiod
1. Thelikelihoodratiotest(LRT)andrelatedF-test.
either.
2. Bayesfactors.
AlldatahavebeenobtainedfromtheUKSwift archive7(Tyler
3. Bayesianposteriorpredictiveprobability.
etal.2006)andprocessedthroughxrtpipelineversion0.10.3,8
4. Monte Carlo test for peaks in data after ‘‘matched filter’’
usingversion008calibrationfilesandcorrectingfortheWTmode
smoothing.
gainoffset(if present).Version008oftheCALDBisamarked
ExamplesoftheapplicationofthesemethodstoGRBX-rayspec- improvementoverthepreviousrelease(Campanaetal.2006a),9
tracanbefoundinFreemanetal.(1999),Yoshidaetal.(1999), however,itmaybethecasethatlow-energycalibrationfeatures
Piro et al. (1999), Protassov et al. (2002), Rutledge & Sako havestillnotbeenoptimallycorrected.Wedidnotapplyanysys-
(2003),Tavecchioetal.(2004),Butleretal.(2005,2007),Sako tematiccorrectionfactortotheerrorsof ourspectrabecausethe
etal.(2005),andreferencestherein. recommendedfactorisverymuchsmallerthanthestatisticaler-
Inalltheapplicationscitedabove,anunderlyingcontinuum rorsinourspectra.Grade0Y2data,usingextractionregionsof
was assumed, usually in the form of an absorbed power law 20;3pixels,forbothsourceandbackgroundregionshavebeen
(e.g., using the Wisconsin absorption model; see Morrison & used.Atcountratesbelow100countss(cid:3)1WTmodedatadoes
McCammon1983).Thedetectionofalinethenamountstoacom- notsufferfrompile-up(Romanoetal.2006);however,someof
parisonof twomodels:M ,thesimple‘‘continuum’’model,and thetimeintervalsconsideredcontainsufficientfluxtocausepile-up
0
M ,themorecomplex‘‘continuum+line’’model.Thestrength, effects.FollowingRomanoetal.(2006)wehaveexcludedcentral
1
location,and width ofthe emissionline may berestricted oral- regionswhennecessaryasdetailedintheirAppendixA,splitting
lowedtobefreeparameters. the20;3pixelregionintotwo10;3pixelregionsplacedeither
AsdiscussedindepthbyProtassovetal.(2002)andFreeman sideofthecentralexclusionregion.
etal.(1999),therearestrongtheoreticalreasonswhytheLRTis Allspectralfittingandsimulationshavebeencarriedoutusing
notsuitableforassessingthesignificanceof emissionslines,de- XSPECversion12.2.1aborhigherwithbackground-subtracted
spiteitspopularityintheliterature.Wewillnotrepeatthoseargu- spectrabinnedto(cid:4)20countsbin(cid:3)1.Thisbinningpermitstheuse
mentshere.Itisthepurposeof thepresentpapertocomparethe ofthe(cid:3)2minimizationasamaximumlikelihoodmethod.Data
relativemeritsoftheremainingthreemethods,intermsoftheircom- foreachGRBbeingconsideredhavebeentimeslicedwiththe
putationalefficiency,robustness,andsensitivitylimitsbyapply- followingcriteriainmind.
ingallthreemethodstoX-rayspectrafromtheSwiftarchive.This
1. Each spectrum must contain 800Y1600 background-
is a particularly rich archive because of the combination of the
subtractedcounts.Thisisacompromisebetweengoodtimereso-
rapidslewresponseoftheSwiftGRBmission(Gehrelsetal.2004)
lutionandspectralquality.
andthepowerfulX-rayTelescope(XRT;Burrowsetal.2005a).
2. If oneormoreflaresarepresentinthedata,whereverpos-
Theremainderof thispaperisorganizedasfollows.Inx2we
sible(notviolatingcondition1),separatespectrawereextracted
providedetailsof thesampleselectioncriteriaandbasicdatare-
fortherisingandfallingsectionsof theflare,sincespectralevo-
duction. The theoretical basis and practical applications of the
lutionislikelytooccuratthistime(Burrowsetal.[2005b],Zhang
threestatisticalmethodsunderinvestigationaredescribedinde-
[2007],andreferencestherein).
tailinx3.Inx4weapplyallof themethodstoPKS0745-19,a
3. If dataareaffectedbypile-upthenthesetimeperiodswere
knownlineemittingsource,todemonstratetheexpectedoutcome
extractedseparatelyfromthenonYpiled-updata.
whenalineispresent,andx5discussesasimulationstudytoas-
sessthelinedetectionlimitsofeachmethodfortypicalSwiftXRT Therangeof800Y1600countswaschosentoensuregoodtime
data.Inx6wediscusstheresultsfromtheSwiftarchivalGRBaf- resolutionwhilemaintainingsufficientcountstoobtainareason-
terglowdata,highlightingseveralGRBswithpotentialadditional ablespectrum,with(cid:1)40Y80spectralbins(eachwith(cid:1)20counts)
spectralcomponents.x7isdedicatedtoadiscussionofourresults over the useful bandpass. The data considered here were taken
andtheircomparisontootherrecentlinedetectionsinthelitera- duringtheearly,brightphasesoftheafterglowevolution(typi-
ture.Finally,x8presentsourconclusions. cally T þ0 s to T þ500 s), during which the X-ray flux and
2. DATA REDUCTION 7 Seehttp://www.swift.ac.uk/swift_live/obscatpage.php.
8 Releasedate2006March16.
Thispaperreportsontheanalysisof windowedtiming(WT) 9 See http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/xrt/SWIFT-
modedatafromGRB050128toGRB060510B,coveringatotal XRT-CALDB-09.pdf.
No. 1, 2008 LINE SEARCHES IN SWIFT X-RAY SPECTRA 589
(possibly)spectrumareoftenhighlyvariable.Timeresolutionis detailedreview).Inthepresentcontextwehavenostrongtheo-
therefore important to reduce the effects of flux/spectral varia- reticalgroundstopreferoneortheothermodel(lineornoline)
tiononthemodelingofindividualspectra.Furthermore,previous andsoassignequalpriorprobabilitiestoourtwomodels.Thus
claimsof emissionslineshaveoftenreportedthefeaturesastran- theratioof thepriorsinequation(2)issettounityandthepos-
sientandsotimeresolutionmaybeimportantfordetectinglines. terioroddsareequaltotheBayesfactor.Inthefollowingweuse
thetermsposteriorodds,odds,andBayesfactorsinterchangeably.
3. ANALYSIS METHODS
The likelihood functions in equation (2) arefunctions of M
i
Asnotedinx1severaldifferentmethodshavebeenusedinthe only.Ifthemodelscontainnofreeparameters(i.e.,arecompletely
pasttoassessthesignificanceoflinedetectionsintheX-rayspec- specified)thenequation(2)canbeuseddirectly.However,if the
traofGRBs.Thethreemethodsthatarethesubjectofthepresent model does contain free parameters, the likelihood will be a
paperarediscussedindividuallyinthefollowingsubsections.The functionof theparametervalues.Inthepresentcontext,where
readerinterestedonlyintheapplicationof thesemethodsmay theparticularvaluesoftheparametersarenotthesubjectofthe
wishtoskiptox4. investigation,theparametersarereferredtoas‘‘nuisance’’pa-
rameters.Inordertoremovethedependenceonthesenuisance
3.1. Bayes Factors
parametersthelikelihoodfunctionmustbewrittenafunctionof
Thegoalofscientificinferenceistodrawconclusionsabout theNparameters(denoteda¼f(cid:4)1;(cid:4)2; : : :;(cid:4)Ng)andintegrated,
the plausibility of some hypothesis or model, M, based on the or ‘‘marginalized,’’ over the prior probability density function
available data D¼fx ;x ;: : :;x g, given the background in- (PDF)fortheparameters:
1 2 N
formationI(suchasthedetectorcalibration,statisticaldistribu-
Z Z
tionof thedata,etc.).However,whenpresentedwithdataitis p(DjM)¼ p(D;ajM)da¼ p(Dja;M)p(ajM)da: ð3Þ
usuallynotpossibletocomputethisdirectly.Whatcanbecalcu-
lateddirectlyinmanycasesisthesamplingdistributionfordata
Themarginallikelihoodisobtainedbyintegratingoverallpa-
assumingthemodeltobetrue,p(DjM;I).Thisisusuallycalled
rameter values the joint PDF for the data and the parameters.
thelikelihoodwhenconsideredasafunctionof MforfixedD.
This jointPDFmaybeseparatedintotheproductof twoterms
Statementsaboutdataconditionalonthemodelmayberelatedto
usingtherulesof probabilitytheory:p(Dja;M)isthelikelihood
statements about the model conditional on the data by Bayes’
theorem.10InitsusualformBayestheoremrelatesthelikelihood functionofthedataasafunctionofthemodelanditsparameters,
andp(ajM)isthepriorforthemodelparameters.Oncetheseare
totheposteriorprobabilityof themodelMconditionalondataD
assignedonecancomputethenecessarylikelihood(afunctionof
(andanyrelevantbackgroundinformationI),writtenp(MjD;I):
Malone)byintegration.TheBayesfactorformodelM (withpa-
1
p(DjM;I)p(MjI) rametersa1)againstmodelM0(withparametersa0)maynowbe
p(MjD;I)¼ : ð1Þ written
p(DjI)
R p(Dja ;M )p(a jM )da
Tdehsectreibrmespo(uMrjkIn)oiwstlheedg‘‘epr(ioorripgrnoobraabnicliet)y’o’fotfhtehemmodoedleplrMioarntod B10¼ R p(Dja10;M10)p(a10jM10)da10: ð4Þ
considerationofthedata(oftencalledsimplytheprior).Theterm
p(DjI) is effectively a normalization term and is known as the Theissuesof howtherelevantlikelihoodsandpriorsareassigned,
priorpredictiveprobability(itdescribestheprobabilitywithwhich andtheintegralscomputed,arediscussedbelow.
onewouldpredictthedatagivenonlypriorinformationaboutthe
3.1.1.Application to High-CountX-RaySpectra
model).Foramoregeneraldiscussionof BayestheoremseeLee
(1989),Loredo(1990,1992),Sivia(1996),Gelmanetal.(1995), Inthelimitof alargenumberof countsperspectralbin,the
andGregory(2005),andfordiscussioninthecontextofGRBline Poisson distribution of counts ineach binwill convergetothe
searchesseeFreemanetal.(1999,theirx3.1.2)andProtassovetal. Gaussiandistribution,andinthiscaseequation(3)canbewritten
(2002).Intherestofthispaperwedroptheexplicitconditioningon intermsof thefamiliar(cid:3)2fitstatistic(Eadieetal.1971),
background information I, but it is taken as accepted that ‘‘no
probabilityjudgmentscanbemadeinavacuum’’(Gelmanetal. L¼ln½p(Dja;M)(cid:5)
1995).
Onesimplewaytorepresenttheposteriorprobabilitiesfortwo ¼(cid:3)1XN ln(cid:1)2(cid:5)(cid:2)2(cid:2)(cid:3)XN ½xi(cid:3)(cid:6)i(a)(cid:5)2
alternativemodelsisintermsoftheirratio,theposteriorodds(see 2 i 2(cid:2)2
i¼1 i¼1 i
Gregory2005,theirx3.5).Thiseliminatesthep(D)term(which
¼constant(cid:3)(cid:3)2=2; ð5Þ
hasnodependenceonM).Ifwedefinetwocompetingmodels,
suchasonewithaline(M )andonewithout(M ),wemaycom-
1 0
putetheposteriorodds: where(cid:2)i istheerroronthedata(e.g.,counts)measuredinthe
ithchannel,and(cid:6)(a) isthepredicted(e.g.,modelcounts)inthe
i
p(M jD) p(M )p(DjM ) p(M ) channelbasedonthemodelwithparametervaluesa.Thelast
O10¼ p(M10jD) ¼ p(M10)p(DjM10) ¼ p(M10)B10: ð2Þ equalitycanbemadesincethetermPln(2(cid:5)(cid:2)i2)isconstantgiven
datawitherrors(cid:2).Thisiswhy,inthehighcountlimit,findingthe
i
Highoddsindicategoodevidencefortheexistenceofalinein parametervaluesatwhich(cid:3)2isminimizedisequivalenttofinding
thespectrum.Thefirsttermontherighthandsideistheratioof themaximumlikelihoodestimates(MLEs)of theparameters,aˆ.
thepriors,thesecondtermistheratioofthelikelihoodsandis
3.1.2.Approximatingthe Posterior
oftencalledtheBayesfactor(seeKass&Raftery[1995]fora
Ingeneraltheintegralsof equations(3)and(4)mustbecom-
10 ForgeneralreferencesrelatingtoBayesiananalysisseehttp://www.astro putednumerically,usingforexampleMarkovchainMonteCarlo
.cornell.edu/staff/loredo/bayes. (MCMC)methods(Gelmanetal.1995,chapter11;Gregory2005,
590 HURKETT ETAL. Vol. 679
chapter 12), which is computationally demanding. However, in van Dyk et al. (2001). The MCMC method does not use an
maximumlikelihoodtheorysaysthattheMLEwillbecomemore analyticalapproximationfortheposterior,andthereforeisamore
Gaussianandof smallervarianceasthesamplesize(numberof generalmethod,butiscomputationallydemanding.Figure1illus-
counts)increases,evenif themodelisnonlinear(chapter11of trates the two posterior distributions calculated for the specific
Gregory2005).Therefore,withsufficientcountsthelikelihood caseof aspectrumfromGRB060124.TheGaussiandatawere
functionwillapproachamultidimensionalGaussianwithapeak computed from 105 random draws from a multidimensional
attheMLElocationaˆ,i.e.,thelocationof thebestfit(minimum Gaussianwith acovariancematrixevaluatedasthe minimum
(cid:3)2)inparameterspace.Furthermore,if thepriorfunctionisrela- (cid:3)2locationusingXSPEC.Thenon-Gaussiandataweregenerated
tivelyflataroundthepeakoftheGaussianlikelihoodwemayap- from105drawsgenerated12bytheMCMCroutineof vanDyk
proximatethepriorterminequation(3)byaconstant,namelyits etal.(2001).Itisclearthatthetwodistributionsarenotidentical
valueatthebest-fitlocation,p(aˆjM).Puttingthistogethermeans butareverysimilarbothintermsofsizeandshape.Inthepresent
wemayapproximatetheposteriorasaGaussian—oftencalled contextitisimportantthatthe‘‘credibleregions’’occupysimilar
theLaplaceapproximation—whichgreatlysimplifiestheintegra- volumesof parameterspace.
tionsinequations(3)and(4),sinceamultidimensionalGaussian The above analyses demonstrate that the Gaussian approxi-
maybeevaluatedanalytically,onceitspeaklocationsandco- mationisreasonablefortheposteriorof thesimple‘‘continuum’’
variancematrixareknown,whichavoidstheneedforcompu- modelM , which isthe denominator ofequation(4).Thesame
0
tationallyexpensivenumericalintegration. willbetrueof themorecomplex‘‘continuum+line’’modelM
1
Theintegralof anun-normalizedmultidimensionalGaussian whenthelineiswelldetected(seex11.3ofGregory2005).When
is(2(cid:5))N/2½det(s2)(cid:5)1/2 timesthepeakvalue,wheres2 istheco- thelineisweaklydetectedtheposteriorwillbeclosetothebound-
variancematrix11evaluatedatthepeak(best-fitlocation)andN aryoftheparameterspace,inwhichcasetheGaussianapproxi-
isthenumberof parameters.Wecannowrewriteequation(3)as mationwillnotbesoaccurate.Indeed,whentheMLEof theline
normalizationisclosetotheboundarythelikelihood(andthere-
pðDjMÞ¼exp(cid:3)(cid:3)(cid:3)2 =2(cid:4)ð2(cid:5)ÞNi=2qffidffiffieffiffitffiðffiffisffiffiffiffi2ffiffiÞffiffip(cid:1)aˆ jM(cid:2); ð6Þ foreposterior)enclosedintheallowedregionofparameterspace
i (i) i i i will be smaller than that given by the Laplace approximation,
whichassumesthattheGaussianfunctionextendstoinfinityin
whichinvolvesthepriordensityonlyatthemodeof thep(aˆ jM) alldirections.Thiswillalsohappen,forexample,whenthebest-
i i
likelihood(i.e.,thedensityattheMLEposition)foreachmodel. fittinglineenergyisnearthelimitof theallowedenergyrange.
Thiscanbesubstitutedintoequation(2)togivetheBayesfactor InsuchcasestherewillbeatendencytooverestimatetheBayes
(seeGregory2005,chapters10Y11) factor (i.e.,favor M ). Butwhen thelineis weakthere may be
1
multiple peaks in the likelihood (and posterior) which are not
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
B ¼exp(cid:1)(cid:3)(cid:1)(cid:3)2=2(cid:2)ð2(cid:5)Þ(cid:1)N=2 det(cid:1)s12(cid:2) p(cid:1)aˆ1jM1(cid:2); ð7Þ afocrceoutrnetaetdthfoerceaxlpculilcaittelydiBnathyeesLfaapcltaocresaopnplryoxasimaartoioung.hWgeutihdeerteo-
10 qffidffiffieffiffitffi(cid:1)ffiffisffiffiffiffi2ffiffi(cid:2)ffiffi p(cid:1)aˆ0jM0(cid:2) thepresenceof aspectralline.
0
3.1.4.AssigningPriors
where(cid:1)N ¼N (cid:3)N and(cid:1)(cid:3)2 ¼(cid:3)2(cid:3)(cid:3)2.Thiscanbecalcu-
lated,usingthea1pprop0riatevaluesof1(cid:3)2and0thecovariancema- Bayesfactorsaresensitivetothechoiceof priordensity.As
trixs2evaluatedatthebest-fitlocationforeachmodel,oncewe statedabove,usingtheLaplaceapproximationtheresultingBayes
assign prior densities to each parameter. See Gregory (2005), factorsaresensitivetothepriordensitiesonlyattheMLEparam-
chapters10Y11,foradiscussionofessentiallythesamemethod. etervalues,butwemustexercisecareinassigningpriordensity
functionsinorderthatthesevaluesarereasonable.Fortunately,the
3.1.3. Validity of theLaplace Approximation priordensitiesforallparametersthatarecommontoM andM
0 1
(suchasphotonindexandnormalization)arethesameforM and
Thereareanumberof waystocheckthevalidityof thisas- 0
sumption.Oneistoinspecttheshapeofthe(cid:3)2(a)surface,which M1,andthereforecanceloutintheratio.Fortheotherparameters
is related to the likelihood surface by L(cid:1)(cid:3)(cid:3)2/2. A Gaussian wehavenocogentinformationexceptfortheirallowedranges.In
likelihoodisequivalenttoaparaboloidallog-likelihood,or(cid:3)2(a). suchcasesweshouldusethe‘‘leastinformative’’priordensities
Ifthecontoursof(cid:1)(cid:3)2appearparaboloidalaroundtheminimum, (see,e.g.,Loredo[1990],Sivia[1996],Gregory[2005],andrefer-
encesthereinforfurtherdiscussion).
inoneandtwodimensions,foreachparameterorpairof param-
Therearewiderangesof possiblelineenergiesandredshifts
eters, the likelihood surface must be approximately Gaussian
andsothelineenergy,E isonlyconstrainedtoliewithinthe
(theconditionalandmarginaldistributionsof aGaussianarealso line
Gaussian,soslicesthroughthe(cid:3)2spaceshouldalsobeparab- useful XRT bandpass, typically 0.3Y10 keV. We therefore as-
signed a uniform prior density p(E jM )¼1/(E (cid:3)E )
oloidal). Thiswasgenerallytrueforthecontinuummodelfor line 1 max min
overthisrange.FormostspectralfitsthelinewidthW wasini-
theXRTdata.
tiallyheldfixedatavaluebelowtheinstrumentalresolution,13
Asafurthertestof theGaussianapproximationwecompared
andlaterallowedasafreeparameter.Forthosemodelsinwhich
the posterior calculated using the Laplace approximation with
thewidthofthelinewasafreeparameter,thewidthwasassigned
the posterior calculated using the MCMC algorithm discussed
auniformpriorovertheallowedrange(usually0.0Y0.7keV):
p(WjM )¼1/(W (cid:3)W ).
1 max min
11 Thecovariancematrixisthesquare,symmetricmatrixcomprisingtheco-
variancesofparameters(cid:4)iand(cid:4)jaselement(cid:2)i2j.Bysymmetry(cid:2)i2j ¼(cid:2)j2i.Thedi-
agonalelementsarethevariancesoftheparameters.Thecovariancematrixmay 12 FollowingvanDyketal.(2001),wegeneratedfiveseparatechains,start-
beestimatedasminustheinversetheHessianmatrixlistingallthesecondde- ingfromdifferent,‘‘overdispersed’’positionswithintheparameterspace(allout-
rivativesoftheloglikelihoodfunction(99L)ij¼@2L/@(cid:4)i@(cid:4)j.GiventhatL¼ sidethe99%confidenceregioncalculatedusing(cid:1)(cid:3)2)andusedtheRˆ1/2statisticto
log½p(Dja;M)(cid:5)¼constant(cid:3)(cid:3)2/2(inthelimitofmanycountsperchannel)the assesstheirconvergence.WecollecteddatafromthechainsonlyafterRˆ1/2<1:01.
covariancematrixmaybeestimatedusing½(cid:2)2(cid:5) ¼2½(99(cid:3)2Þ(cid:3)1(cid:5) .Thesecond 13 (cid:2)¼59 eV (at 5.895 keV) at launch (A. Beardmore 2007, private
ij ij
derivativesofthe(cid:3)2(a)functioncanbeevaluatednumerically. communication).
No. 1, 2008 LINE SEARCHES IN SWIFT X-RAY SPECTRA 591
Inordertotestthedependenceof theresultstothepriorden-
sities,twononinformativepriorassignmentsweremadefortheline
strength(normalization),A,followingthediscussioninGregory
(2005).First,followingx4.2ofSivia(1996)thelinestrengthwas
assignedauniformpriorbetweenzeroandsomeupperlimitA .
max
Previousreportsofemissionlineshaveestimatedthelineflux A
tobeaslittleasafewpercent(Reevesetal.2002;Watsonetal.
2003)orasmuchas(cid:1)40%Y80%(Yoshidaetal.1999;Piroetal.
2000)ofthetotalflux.WeconservativelytakeA tobethetotal
max
flux of the spectrumoverthe evaluated bandpass(i.e.,our con-
straintisthatthelinefluxisbetween0%and100%of thesource
flux).However,therearestrongarguments(Loredo1990;Gelman
etal.1995;Gregory2005)thatsucha‘‘scale’’parametershouldbe
givena‘‘Jeffreysprior,’’p(AjM )(cid:1)1/A,whichcorrespondstoa
1
constantdensityinlog(A).Formallythisisanimproperprior(can-
notbenormalizedsuchthatitsintegralisunity),butonecanapply
reasonableupperandlowerboundsinordertoformaproperprior
density.Followingequation(3.38)of Gregory(2005),weused
p(AjM )¼1/Aln(A /A ).InthepresentcontextA /A ¼
1 max min max min
800,sinceareasonablelowerlimittotheX-raycountsfromaline
isonecount,andareasonableupperlimitis800,thetotalnumber
of counts in the spectrum. This yields p(AjM )¼1/6:68A as a
1
normalized Jeffreys prior. The prior density is therefore higher
forweakerlinesintheJeffreyscasecomparedtotheuniformcase
atvalues(i.e.,over1:25;10(cid:3)3A (cid:6)A<0:15A ),andlower
max max
forstrongerlines.
Theratioofthepriordensitiesatthemodesof thetwolikeli-
hoodfunctionsisthensimply
p(aˆ jM )
1 1 ¼p(Eˆ jM )p(Wˆ jM )p(AˆjM )
p(aˆ jM ) line 1 1 1
0 0
1
¼ ð8Þ
½ðE (cid:3)E ÞðW (cid:3)W ÞA (cid:5)
max min max min P
in the ranges E 2½E ;E (cid:5), W 2½W ;W (cid:5), and zero
line min max min max
elsewhere.Here,A ¼A intheuniformcaseorA ¼6:68Aˆ
P max P
intheJeffreyscase.WehaveusedbothuniformandJeffreyspriors
intheanalysisdiscussedbelow(seex6.11).
3.2. Posterior Predictive p Values (ppp)
Theuseofposteriorpredictive pvalues(ppp)wasadvocated,
anddemonstratedbyapplicationtoGRBspectra,byProtassov
etal.(2002;seetheirx4.1foradescriptionof theirmethodand
theirx5foritsapplicationtoGRB970508).LikeBayesfactors
thismethodisgroundedinBayesianprobabilitytheory.
Oneusestheposteriordensity, p(ajD),forthemodelparame-
ters that are conditional on the data, which defines our state of
knowledgeabouttheparametersgiventhedataandtheavailable
priorinformation,todeterminetheposteriorpredictivedistribu-
tion, which is the distribution of possible future data predicted
basedontheobserveddata.(‘‘Predictive’’becauseitpredictspos-
siblefuturedatasetsand‘‘posterior’’becausetheparametersare
drawnfromtheposteriordensityoftheparameters.)Theposterior
predictivedistributionis
Fig.1.—Marginalposteriordistributionsforthecontinuumparametersofan Z
absorbedpowerlawfittoanXRTspectrumofGRB060124.Thecontoursen- p(cid:1)DsimjD(cid:2)¼ p(cid:1)Dsim;ajD(cid:2)da
close80%,50%,20%,10%,and5%ofthedistribution(andthereforecorrespond
to20%,50%,80%,90%,95%credibleregionsfortheparameters).Thefilledcon- Z
tourswerecomputedassumingtheposteriorisaGaussian.Theopencontourswere ¼ p(cid:1)Dsimja(cid:2)pðajDÞda; ð9Þ
computedusingMCMCsimulationsfromtheroutineof vanDyketal.(2001).The
twodistributionsareclearlyverysimilar.Seex3.1.3fordetails.
whereDsimrepresentsthepossiblefuturedatasets(simulations).
Inpracticetheposteriordensityisusedtogenerateasetofrandom
parametervaluesasim(i¼1;2; : : :),andeachoftheseisusedto
i
592 HURKETT ETAL. Vol. 679
simulatearandomdataset Dsim.Thesetof simulateddatafrom For thepurposesof thepresentpaper weuse astheteststa-
i
allthepossiblerandomparametersdefinestheposteriorpredictive tisticthechangeinthe(cid:3)2fitstatistic14betweenthetwomodels,
distributionforsimulateddata.Thisinturncanbeusedtodefine M and M . This is equivalent to the formulation discussed in
0 1
the posterior predictive distribution for some test statistic T(D) Protassov et al. (2002). The observed data were fitted with the
(whichisafunctionof thedata): modelM andthecovariancematrixevaluatedatthebest-fitpoint
0
usedtoconstructthemultivariateGaussiandistributionfromwhich
phT(cid:1)Dsim(cid:2)(cid:6)(cid:6)Di¼Z phT(cid:1)Dsim(cid:2)(cid:6)(cid:6)aipðajDÞda ð10Þ parametervalueswererandomlydrawn.15Foreachsetof model
(cid:6) (cid:6) M parametervaluesaspectrumwassimulatedwiththeappro-
0
priate response matrix and exposure time, with counts in each
channeldrawnfromaPoissondistribution,andbinnedinthesame
(comparewitheq.[9]).Thepppisthefractionof thisdistribu-
mannerastheobserveddata.
tionforwhichT(Dsim)>T(D),i.e.,theareaofthetailofthedis-
In order to calculate the test statistic for each simulation,
tribution with valuesoftheteststatisticmore extreme thanthe
T(Dsim),itwasnecessarytofiteachsimulateddatasetwiththe
valuefromtheobserveddata; i
twocompetingmodelsM andM ,findthebest-fittingparame-
0 1
p¼Z 1 phT(cid:1)Dsim(cid:2)(cid:6)(cid:6)DidDsim; ð11Þ tceormsfpourtaetaicohnaolnlye,eaxnpdecnosimvepumteul(cid:1)tid(cid:3)ii2m.eTnhsiisonneaclepsasraarmilyetienrveosltvimesaa-
(cid:6)
T(D) tionforeachof theNsimulations.WeuseasstandardN ¼104
simulations,whichyieldsa pvalueaccuratetofourdecimalplaces
wheretheintegrationistakenovertheposteriorpredictivedis- atveryhighestandlowest pvalues(thereisanuncertaintyonthe
tributionof Dsim.AssuchthepppvalueisaBayesiananalogof pppvaluefromthefinitenumberof simulationswhichisroughly
the pvalueof nullhypothesistestsfamiliarfromclassicalstatis- ½ p(1(cid:3)p)/N(cid:5)1/2 fromthebinomialdistribution).Thisisaccept-
tics(e.g.,theFor(cid:3)2tests).Seechapter6ofGelmanetal.(1995) ablefordeterminingpvaluesaslowas p(cid:1)10(cid:3)4,i.e.,99.99%
orGelmanetal.(1996)forageneraldiscussionofthepppmethod, ‘‘significance.’’
andProtassovetal.(2002)forapplicationtoGRBdata. AsafurthertestofthevalidityoftheGaussianassumptionfor
Using theposteriorpredictivedistribution fromequation(9) theposterior(seealsoxx3.1.2Y3.1.3)wehavecomparedresults
onecanproducealargenumberofrandomsimulateddatasetsto withandwithoutthisassumption.Inparticular,wecalculatedthe
beusedinaMonteCarloschemetocalculatetheintegralofequa- pppvalueforaspectrumofGRB060124usingGaussianparam-
tion(11)numerically.ThestepsforaMonteCarlomethodfor etervaluesandalsousingvaluesgeneratedbytheMCMCmethod
computingtheposteriorpredictivedistributiontocalibratethe discussedbyvanDyketal.(2001).Thetworesultswerereason-
teststatisticT areasfollows. ably close (p¼0:050(cid:7)0:007 from the Gaussian simulations
and p¼0:073(cid:7)0:008fromtheMCMC,basedon103simula-
1. Computethevalueoftheteststatisticfortheobserveddata,
tions).Thisconfirmsthepointmadeinx3.1.3,thattheGaussian
T(D).
assumptionisreasonableforthesedata.
2. Randomly draw N sets of M model parameter values a
0 i
for i¼1;2;: : :;N according to the appropriate posterior dis- 3.2.2. Automated FittingofGRB Spectra
tributionp(ajD).
3. Foreachofi¼1;2;: : :;NsimulateadatasetDsimusing Giventhenumberof simulateddatasetsonemustresorttoan
therandomlydrawnparametervaluesa.Thisaccountisforun- automatedmodelfittingprocedure.Thishasitself beenthecause
i ofsomedebate,withsomeauthors(e.g.,x5of Rutledge&Sako
certaintiesintheparametervalues.
2003)claimingthatautomaticroutinesdonotrobustlyfindthe
4. Foreachofthesimulateddatasetscomputetheteststatistic
T(Dsim). This is the posterior predictive distribution of the test best-fittingparametervalues(minimum(cid:3)2).Thealgorithmused
i byXSPECfor(cid:3)2minimizationistheLevenberg-Marquardtal-
statisticgiventheobserveddataD.
gorithm,whichisefficientandveryeffectivewhenthe(cid:3)2space
5. Computetheposteriorpredictivepvalueasthefractionof
iswellbehaved(e.g.,withonlyonelocalminimum).However,
simulated data sets that gave a test statistic more extreme than
asthisisa‘‘local’’routine,thereisnoguaranteeof findingthe
thatfortheobserveddata,
‘‘global’’minimumin(cid:3)2,anditispossiblethattheresultsarebi-
asedbythepresenceofotherlocalminima.Forthepresentpaper
p¼ 1 XN (cid:2)(cid:7)T(cid:1)Dsim(cid:2)(cid:3)TðDÞ(cid:8); ð12Þ wehaveemployedseveraladditionstothestandardLevenberg-
N i Marquardt minimization algorithm in order to mitigate these
i¼1
problems.
Oncealocalminimumin(cid:3)2isfoundthesurroundingregion
where (cid:2) is the Heaviside step function which simply counts
instanceswhereT(Dsim)>T(D). of parameterspaceisexploredforsignsof otherminima.Each
i parameterinturnhasitsvalueincreasedanddecreaseduntilthe
Thenumberofsimulations,N,mustbelargetoensureagood (cid:3)2 is increased by at least (cid:1)(cid:3)2 ¼2:7, while simultaneously
approximationtotheintegralof equation(11)(whichisamul-
tipleintegral,beingitselftheintegralofthefunctioncomputedby
eq.[10]).SeeProtassovetal.(2002)formoredetaileddiscussion. 14 The(cid:1)(cid:3)2statisticisfamiliartomostX-rayastronomersandwasusedinthe
Bayesfactorsmethodabove.Herewenotethatitisequivalenttothelikelihoodra-
3.2.1. Application toGRB X-RaySpectra tiotest(LRT)statistic,sinceusingeq.(5)wehave(cid:1)(cid:3)2¼(cid:3)2 (cid:3)(cid:3)2 ¼(cid:3)2lnk,
(0) (1)
wherek¼p(Djaˆ0;M0)/p(Djaˆ1;M1)istheratioofthelikelihoodmaximaof the
Asdiscussedabove,wemayapproximatetheposteriordensity twomodels.UndertheassumptionsforwhichtheLRTisvalidthisshouldbedis-
fortheparametersusingamultidimensionalGaussiancenteredon tributedas(cid:3)2withdegreesoffreedomequaltothenumberofadditionalfreeparam-
theMLEvaluesandwithashapedefinedbythecovariancematrix etersinmodelM1comparedtoM0.ThereasonforchoosingtheLRToverrelated
evaluated at the peak (s2). The randomized parameter values statistics,suchastheF-test,isthatLRTismorepowerful.SeeFreemanetal.(1999),
Protassovetal.(2002),andreferencesthereinfordetails.
neededforstep2abovemaythenbegeneratedwiththeCholesky 15 Inpracticethiswasperformedusingthetcloutsimparscommandin
method. XSPEC.
No. 1, 2008 LINE SEARCHES IN SWIFT X-RAY SPECTRA 593
allowingtheotherparameterstovaryinordertominimize(cid:1)(cid:3)2. 3.4.Comparison of the Methods
If anynon-monotonicityin(cid:3)2isdetectedduringthissearchthe
Thethreemethodsdiscussedabovehavedifferenttheoretical
volumeof parameterspaceexploredisincreasedbyincreasing
motivations and underlying assumptions and require different
the value of (cid:1)(cid:3)2. If during the course of this search (cid:1)(cid:3)2 be-
amountsofcomputingpower.TheBayesfactormethodisbased
comesnegative(meaningthereisalowerminimumnearby)the
on a simple application of Bayes theorem combined with the
Levenberg-Marquardtalgorithmisrestartedfromthepositionof
Laplaceapproximationandassumesuniformpriorsonthemodel
thisnewminimum.Theentireprocessisrepeatedbyperturbing
parameters (or Jeffreys prior for the linenormalization).Asdis-
eachparameterinthiswayuntilnofurtherimprovementcanbe
cussedabove,thismaynotbetheoptimalassignment.However,
madebytheadjustmentof anyof them.
despite its possible drawbacks, the simple priors and Laplace
Theabsorbedpowerlawmodel(M )hasonlythreeparameters
0 approximationmakethecalculationextremelysimple,requiring
(photonindex(cid:3),normalization,andabsorptioncolumndensity),
onlytheevaluationofequations(7)and(8),whichrequiretheval-
andinallcasesfindingthe(cid:3)2minimumwasstraightforwardus-
uesof(cid:3)2andthecovariancematricesforthebest-fittinglineand
ingtheaboveprocedure.ThealternativemodelM ,whichin-
1 line-freemodels,anddetailsofthefreeparametersandtheiral-
cludestheemissionline,requiredmorecarebecausethepresence
lowed ranges. As such, it is useful as a ‘‘quick and easy’’ test.
ofalinewithunknownenergymaycauselocal(cid:3)2minimaatdif-
Thedependenceof choiceofpriorsmaybeassessedbycompar-
ferentenergieswithinthewidebandpass.Aninitial‘‘bestguess’’
ingtheresultscomputedusingtheuniformand Jeffreysprior.
lineenergywascomputedforeachspectruminthefollowingway.
Bycontrast,theRSandpppmethodsrequirealargenumber
Anabsorbedpowerlawplusemissionlinemodelwasconstructed
of randomdatasetstobesimulatedandanalyzed,andarethere-
usingthebest-fittingparametersof modelM andaddinganun-
0 foreconsiderablymorecostlyintermsofcomputingtime.There
resolvedemissionlinefixedatsometrialenergyE andvarying
i isnocompellingtheoreticalreasonforapplyingamatchedfilter,
theotherparameters(includingthelinenormalization)tofindthe
asintheRSmethod,althoughitshouldbenotedthatthemethod,
minimum(cid:3)2.OnehundredvaluesofthetrialenergyE wereused,
i asimplementedabove,iscalibratedusingtheappropriateposte-
evenlyspreadovertheentireusefulbandpass,andthevaluethat
riorpredictivedistribution.TheadvantageoftheRSmethodisthat
recordedthelowest(cid:3)2wasselectedasthe‘‘bestguess’’fortheline
nomodelfittingisrequired,whichisoftenatime-consumingpro-
energy.TheenhancedLevenberg-Marquardtalgorithmdescribed
cessandcanleadtobiasedresultsifnothandledproperly(x3.2.2).
above was then used to find the global (cid:3)2 minimum starting
Thepppmethodisgroundedinthetheoryof Bayesianmodel
fromthisposition.Simulationtestsandcomparisonwithinter-
checking(Gelmanetal.1995;Protassovetal.2002)butrequires
active fittingdemonstratedtheautomaticproceduredescribed
time-consumingfitstobeperformedoneachsimulatedspectrum,
abovewasanefficientandveryrobustprocedureforfindingthe
and is therefore the most computationally demanding method
globalminimum.
byaclearmargin.However,itisarguablythemostrigorousin
the sense that it is less sensitive to the choice of priors than
3.3. Rutledge and Sako Smoothing (RS)
areBayesfactors(Gelmanetal.1996;Protassovetal.2002),
Rutledge&Sako(2003)proposedanalternativemethodfor anddoesnotapplyanadhocsmoothing,asintheRSmethod,
linedetectionusinga‘‘matchedfilter’’tosmooththeobserved thatmayactuallyacttosuppressrealspectralfeaturesinsome
countspectrumwiththeaimofremovinglow-significancenoise cases.
andemphasizinganyspectralfeatures.Thedistributionof peak ThesimulationsusedforbothRSandpppmethodsweregen-
fluxesinthesmoothedspectrumisthencomparedtotheresultof erated assuming a Gaussian posterior for the three parameters
MonteCarlosimulationstocalibratetheirsignificance(p value). of M0,which,asdiscussedinx3.1.3,isagoodapproximation.
ThecountsperPHAchannelareextractedfromtheobserved Again,thisapproximationwasmadetoincreasecomputational
X-rayspectrumandthensmoothedusinganenergy-dependent efficiency,sinceGaussiandeviatesaretrivialtogeneratewiththe
kernel(aGaussianhavingaFWHMequaltothespectralresolu- Choleskymethod.InsituationswheretheGaussianapproxima-
tionofthedetector;seeeq.[2]ofRutledge&Sako2003)topro- tionisnotvalidand/orthenumberofspectraissmallenough
ducethesmoothedspectrumC(E).ThedistributionofC(E)is thatconsiderablymorecomputingtimemaybespentoneach,
thencalibratedusingMonteCarlosimulationsofspectragener- thepppmethodorBayesfactorsmaybecomputedusingresults
atedusingthemethoddiscussedinx3.2thatemploysposterior fromMCMCsimulations(vanDyketal.2001;Protassovetal.
predictivedatasets.(WenotethatRutledge&Sako[2003]and 2002),whichallowsforamoreaccurateevaluationof theposte-
Sako et al. [2005] did not randomize the parameter values but riordensity.
usedfixedMLEvaluestogeneratealltheirsimulations.Thisis
3.4.1.AlternativeApproximateMethods
equivalenttoassumingtheposteriortobeadeltafunctionlocated
atthebest-fitpoint,whichisclearlyabadapproximationinmany Thestatisticsliteraturecontainsmanymethodsdevelopedfor
cases.)Eachsimulationisinturnsmoothedusingthesameenergy thepurposeofmodelselection.Inx1welistedfourmethodsthat
kerneltoproduceC(E) .TheC(E) valuesarethensorted havepreviouslybeenappliedtotheproblemoflinedetectionin
sim;i sim
indescendingorderforeachPHAchannelseparately.Thusthe X-raydatafromGRBs.Onemethodthathasnot,toourknowledge,
99thpercentilelimitof theC(E) isthenfoundbyextract- beenappliedspecificallytoGRBlinedetectionistheBayesian
sim;global
ingthe100thhighestvalueof C(E) ineachPHAchannel. information criterion (BIC; Schwartz 1978). This aims to ap-
sim
The smoothed observed spectrum, C(E), is then plotted proximatethelogarithmoftheintegratedposteriorprobability
alongsidethenthpercentilelimits,whichwehavechosenforthis foramodelwithkparametersgivendatawithasamplesizeN.
analysistobe90.00%,99.00%,99.90%,and99.99%.Wherever TheBICtakestheformofthelogarithmofthelikelihoodwitha
C(E)exceedsagivenlimitthenwehavedetecteda‘‘feature’’at penaltyterm:
that confidence limit. Thus a line would show up as a narrow
excesswhileotherthermalemissioncomponentswillshowupas
broadexcess,bothofwhichareeasilydistinguishable. BIC¼(cid:3)ln½pðDja;MÞ(cid:5)þ(k=2)lnN: ð13Þ
594 HURKETT ETAL. Vol. 679
Fig.2.—SpectralfittoPKS0745-19withanabsorbedpowerlawplusanar-
rowGaussianemissionlinemodel.Theredshiftedironlineat6.07keVisclearly Fig. 3.—PKS0745-19:RSmethod.Confidencecontoursmarkthesignifi-
visible.Notealsotheresidualsat0.6keVand2.3keV,whicharethoughttobe canceofthespectralfeatures.Insetsfocusonenergyrangesofinterest.[Seethe
duetotheoxygenandgoldedgesrespectively.[Seetheelectroniceditionofthe electroniceditionoftheJournalforacolorversionofthisfigure.]
Journalforacolorversionofthisfigure.]
slightlyimprovedthefitwith(cid:3)2/(cid:7) ¼610/509.Allspectralpa-
rametererrorsarequotedat90%confidence.
ThemodelwiththesmallestBICisfavored.Thedifference
AddingaGaussiancomponent toanabsorbedpowerlawfit
betweentheBICvaluesfortwocompetingmodels(oftencalled
naturallyproducedasignificantlyimprovedfittothedata((cid:3)2/(cid:7) ¼
the Schwartz criterion) is therefore S ¼(cid:3)lnkþ((cid:1)k/2)lnN
539/508)with(cid:3)¼2:36þ0:04andalineat6:07þ0:02keV(width¼
(seefootnote14),andisaroughapproximationtothelogarithm (cid:3)0:03 (cid:3)0:02
0:06þ0:02keV).ThespectralfittothismodelcanbeseeninFig-
of theBayesfactor(x4.1.3of Kass&Raftery1995). ure2(cid:3).0T:0h3isissupportedbytheBayesfactorof 7;1014forasin-
In the high count (large sample size) limit (see eq. [5]) the
Schwartz criterion becomes 2S ¼(cid:3)(cid:1)(cid:3)2þ(cid:1)klnN. Whether glelinebeingpresent.RSanalysisofthespectrum,Figure3,also
clearlyshowedthepresenceofaGaussianfeatureat(cid:1)6.07keV
ornottheBICformodelM issmallerthanthatforM isthen
1 0
equivalenttothecriterion(cid:1)(cid:3)2 >(cid:3)(cid:1)klnN.Inthepresentcase withasignificancefarinexcessofthe99.99%confidencelimit.
Thepppanalysisplacedasignificanceof >99.99%onthisfeature.
thedataareselectedwithfixedN,and(cid:1)k ¼(cid:3)2fortheaddition
Aninterestingpointtonoteisthatthereareshallow‘‘excesses’’
of afixedwidthline,suchthattheBICisequivalenttoapplying
thesame(cid:1)(cid:3)2criteriontoeachspectrum,mechanicallythesame at(cid:1)0.6and(cid:1)2.3keVintheRSplot(Fig.3)thatareclearlynot
linefeatures.Coherent,low-level,positiveexcessesarealsoseen
astheLRT,althoughwithadifferent(generallyhigher)thresh-
inthespectralfitattheseenergies(Fig.2).Eitherthepowerlaw
oldvalue.Therefore,inthepresentcontexttheapplicationofthe
componentisnotmodelingthedataadequatelyatthesepoints,the
BICwouldbeequivalenttoaslightlymoreconservativeapplica-
energyscaleforthisspectrumhasanoffsetorthecalibrationfiles
tiontheLRT(seefootnote21).However,asnotedinProtassov
are lessaccuratearound thesetwo energies. Applyingthe gain
etal.(2002)andelsewhere,theBICisoftenapoorapproxima-
fitfunctioninXSPECimprovesthemodelfitssignificantlyby
tiontotheintegratedposteriorprobability,andasdiscussedby
addinganoffset16of(cid:3)0.07keV(nochangetotheslope).The
Kass&Raftery(1995)isgenerallyaworseapproximationthan
absorbedpowerlawmodelimprovesfrom(cid:3)2/(cid:7) ¼635/511to
theLaplaceapproximationemployedtocalculateBayesFactors
606/509andthemekalcomponentmodelimprovesfrom(cid:3)2/(cid:7) ¼
inx3.1.2.
610/509 to 585/507. As a result the two shallow ‘‘excesses’’ at
(cid:1)0.6and(cid:1)2.3keVbecomefarlessprominent.
4.RESULTS FROM AN IRON LINE-EMITTING SOURCE
The feature at (cid:1)0.6 keV could be attributed to the detector
Asafirstdemonstrationoftheabovemethodsweappliedthem oxygen absorption edge at 0.54 keV. Applying the (cid:3)0.07 keV
toanon-GRBSwiftdataset.Ideallywewouldprefertoexamine offsetbringsthe(cid:1)0.6keV‘‘line’’inconjunctionwiththisedge,
asourcewithaGRB-likespectrum,withasimilarcountrate,but thusreducingitssignificancebelowthepointatwhichwewould
containingaclearlyidentifiedemissionlinefeature.However, considerittobearealdetection.Wenotethatthe(cid:1)2.3keVfea-
itisdifficulttofindasourcethatmeetsallof thesecriteria.We tureiscoincidentinenergywiththegoldedgeduetotheXRTmir-
chosethePCmodecalibrationdataset(combiningallavailable rors.Wehaveconfirmedthatthisfeature isnotduetoanybad
datafrom2005May10to2005September2)of PKS0745-19 pixelorhotcolumnissues.17
(DeGrandi&Molendi1999;Chen etal.2003).This testhas
5. TESTING THE THREE METHODS
somelimitationsasPKS0745-19isfainterthantheGRBsana-
AND DETERMINING DETECTION LIMITS
lyzedinthispaperanditisobservedinadifferentmode.
PKS0745-19isagalaxyclusterwithathermalspectrumand In this section we discuss the sensitivity limits of the three
aknownlineat6.07keVinSwift’sobservations,whichisared- methods,i.e.,theweakestlinesthatcanbereliablydetectedwith
shifted6.7keVironline(z¼0:1).Eventhoughtheunderlying eachofthethreemethods,forobservationsofthetypediscussed
spectrumisthermal,withmultipletemperaturecomponents,itcan
bemodeledbyanabsorbedpowerlawcontinuumwherethepower
16 Seehttp://swift.gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/xrt/xrt_bias.pdf.
lawindex,(cid:3),is2:34þ(cid:3)00::0033((cid:3)2/(cid:7) ¼635/511).Anadditionalmekal 17 Seehttp://swift.gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/xrt/SWIFT-
(Meweetal.1985;Arnaud1996)component,withkT ¼0:19þ0:03, XRT-CALDB-01_v5.pdf.
(cid:3)0:01
No. 1, 2008 LINE SEARCHES IN SWIFT X-RAY SPECTRA 595
Fig.4.—NarrowGaussianline(width<instrumentalresolution)foraspectrumcontaining800counts(left)and1600counts(right).Comparisonof thedetection
limits,inequivalentwidth(keV),ofthethreemethodsovertheenergybandpassofSwift.Thedataareasfollows.Dottedlines:Bayesfactoranalysis:solidlines:RSmethod;
dashedlines:posteriorpredictivepvalueanalysis.[SeetheelectroniceditionoftheJournalforacolorversionofthisfigure.]
inx2,ofa‘‘typical’’Swifteraburst.Thisisdonebysimulating p(M jD)<0:01wastakenasthecriterionfordetection.Thisis
0
XRTdatawithacontinuumspectralmodeltypicalof theGRBs equivalenttop(M jD)>0:99,andapproximatelyequivalentto
1
observedwithSwift,butincludinganemissionline,andthenap- aBayesfactorB >100,whichisconventionallytakenasstrong
10
plyingthethreemethodsdescribedaboveforlinedetection. evidenceinfavorof M overM (Kass&Raftery1995).How-
1 0
Inordertogeneratethesimulateddataweuseafiducialspec- ever,westressthattheinterpretationofpvaluesandBayesfactors
tralmodelcomprisingapowerlawwithphotonindex(cid:3)¼2:0, arefundamentallydifferent.Apvalueisthetailareaoftheprob-
normalization(at1keV)of N ¼0:9photonskeV(cid:3)1s(cid:3)1andan abilitydensityfunctionoftheteststatistic,assuminganullhy-
absorptioncolumndensityof N ¼1:8;1021cm(cid:3)2(seeTable pothesis(M )istrue,andisusedtodecidewhetherornottoreject
H 0
2of Campanaetal.2006c).Theseparametersaretypicalofthe thehypothesis.Assuch,a pvalueisnottheprobabilityforthe
X-rayspectraof Swifterabursts.18Inordertomeasurethesen- model M ; instead it corresponds to the frequency of more ex-
0
sitivityofthethreedetectionmethodstolinesinXRTdata,spec- tremeteststatistics(e.g.,(cid:1)(cid:3)2)givenalargenumberof repeatex-
traldataweresimulatedusingtheabovemodelplusoneGaussian periments(assumingthenullhypothesis).Bycontrast,p(M jD)
0
emissionline,andsubjectedtoeachof thethreeprocedures.A istheposteriorprobabilityformodelM basedondataDandthe
0
rangeofvaluesforlineenergy,normalization,andintrinsicwidth priors(inthepresentcaseweusedanapproximationthereof),as
wereusedinordertocalibratethedependenceof themethodsto p(M jD)isforM ,andBayesfactorsareusedtoselectbetween
1 1
thelineparameters.19 twomodelsbasedontheratioofthesetwo.Thisfundamentaldif-
The ppp and RS method result in p values with the conven- ferenceintheinterpretationof Bayesfactorsmeansthereisno
tionalfrequentistinterpretation.If wesetthedetectionthreshold expectationthat(cid:1)isthefrequencyof typeIerrorsfromalarge
at (cid:1), and identify a detection as p(cid:6)(cid:1) then the rate of type I number of repeat observations when using a p(M jD)<(cid:1)
0
errors(i.e.,falsepositivedetections)willbe(cid:1).Forthepurpose criterion.
ofsensitivityanalysisweused(cid:1)¼0:01,equivalenttoa‘‘99% Inx3.1.3weconfirmedthatusingtheLaplaceapproximation
significance’’criterion. Incontrast tothese,theBayesfactoris assumptioninthecalculationof theBayesfactorwasvalidfor
theratioofthemarginallikelihoodsofmodelsM andM ;inthe thefiducialabsorbedpowerlawspectralmodel.Thesamewas
1 0
caseofuniformpriorsforthetwomodelsthisistheratioof pos- alsofoundtobetrueofthespectrawithsimulatedGaussianlines
teriorprobabilitiesB ¼ p(M jD)/p(M jD)wheretheprobabili- atandabovethedetectionlimitdetailedabove.
10 1 0
tiesareinterpreteddirectlyasprobabilitiesformodelsM andM , ForeachvalueofthelinenormalizationwecalculatedtheBayes
0 1
respectively. factorsfor50independentsimulationsandcalculatedthep(M jD)
0
Forthepurposeof numericalcomparisonwiththe pvalues, values for each. We then averaged the p(M jD) values at each
0
the Bayes factors were converted into probabilities [assuming normalizationandlinearlyinterpolatedbetweenpointsatadjacent
p(M jD)þp(M jD)¼1;seeeq.(3.19)of Gregory2005],and normalizationvaluestomap p(M jD)asafunctionofnormaliza-
1 0 0
tion.Thelimitingsensitivitywastakentobethenormalizationat
18 Wehaveassumedaredshiftz¼0forthefiducialburstspectrum.The whichthemean p(M0jD)valuefallsbelow0.01.
averageofthemeasuredredshiftsforSwiftGRBsishigherthanthis(seehttp:// Figure4showsthedetectionlimitsforanintrinsicallynarrow
www.astro.ku.dk/~pallja/GRBsample.htmlfortheupdatedvalues).However,it line (W ¼0) at different energies for spectra with (cid:1)800 and
shouldbenotedthatincreasingzcausestheeffectsofabsorptionbythehostgal- (cid:1)1600counts(leftandrightpanels,respectively).Thelimiting
axyabsorption(whichtendstodominatethetotalabsorptioncolumn)toshiftout
sensitivitiesareshowninunitsof equivalentwidth(keV),which
oftheobservedbandpass,meaningthereisrelativelymorefluxatlowerenergies
(<1keV).Thecalculateddetectionlimitsshouldberepresentativeof Swiftera iseasiertointerpretphysicallythantheabsolutenormalization,by
burstsalthoughperhapsconservativeatlowerenergies. comparingthenormalizationtotheunderlyingcontinuummodel.
19 Therangesofvaluesusedforthelinesimulationsareasfollows.Normali- Figures5and6showthedetectionlimitsfordifferentlinewidths
zationsof1;10(cid:3)7!100photonscm(cid:3)2s(cid:3)1takeninlogarithmicallyincreasing
(W ¼0:2and0.7keV,respectively).TheBayesfactorpointsin
steps,lineenergiesof 0.4,0.6,0.8,1.0,2.0,3.0,4.0,5.0,7.0,and9.0keV,and
thesefigureshavebeencalculatedusingtheuniformprior,rather
intrinsicwidthsof 0.0keV(i.e.,unresolved),0.2keV(broadline)and0.7keV
(broadcontinuumexcess). thanthe Jeffreysprior.Seex6.11forfurtherdiscussiononthe
596 HURKETT ETAL. Vol. 679
Fig.5.—BroadGaussianline(width=0.2keV)foraspectrumcontaining800counts(left)and1600counts(right).Comparisonofthedetectionlimits,inequivalent
width(keV),ofthethreemethodsovertheenergybandpassofSwift.Thedataareasfollows.Dottedlines:Bayesiananalysis;solidlines:RSmethod;dashedlines:ppp.
[SeetheelectroniceditionoftheJournalforacolorversionofthisfigure.]
effectof usingthetwodifferentpriorsinthecalculationof the of0.01.Similarly,theselibrarieswereusedtocomputethe99.00%
Bayesfactorsfortheobserveddatasets. significancecontourfromthefiducialmodelfortheRSmethod.We
ThepppandRSmethodsrequirealargenumberof spectral pointoutherethatthesesimulationlibrarieswereusedonlyforthe
simulationsinordertocalibratetheirdistributionandestimatethe purposesofcomparingthedifferentalgorithms.Fortheanalysisof
pvalue.Thecomputationaldemandsofthis20aresuchthatitwas realobservations(discussedbelow),eachobservationwasassessed
notfeasibletoproduceasufficientlylargesetofsimulationsto usingindependentlygeneratedsimulationsmatchingtheparticular
carryoutthemethodsoneachandeverylinespectrum(e.g.,which observationalparameters.
includesseveralspectraateachtrialvalueof lineenergy,width,and Forthepppmethodeachof thespectracontainingalinewas
normalization,forboth(cid:1)800and(cid:1)1600countspectra).Wethere- fittedwithanabsorbedpowerlawwithandwithoutanadditional
foreconstructedtwolibrariesof104simulations,onefor(cid:1)800and Gaussiancomponent,andthechangein(cid:3)2noted.The(cid:1)(cid:3)2val-
onefor(cid:1)1600countspectra,thatcouldbeusedforeachtest.These ueswereaveragedateachnormalization,andthesepointslinearly
wereconstructedbysimulatinganappropriatespectrumbasedon interpolated,tomapthe(cid:1)(cid:3)2asafunctionof normalization.As
thefiducialmodel,andusingthistogeneratetheposteriorpredic- withtheBayesfactor,thelimitingsensitivitywastakentobethe
tivedistributionfromwhichtodraw104simulationsfollowingthe normalizationatwhichthemean pvaluefallsbelow0.01,calcu-
recipediscussedinx3.1.2.Theselibrarieswerethenusedtocali- latedusingtheappropriatedvalueof(cid:1)(cid:3)2valuefromeachsimu-
bratethedistributionof the(cid:1)(cid:3)2statisticforthepppmethodand lation library. The limitingsensitivity as a function ofenergy is
thustocalculatethevalueof (cid:1)(cid:3)2thatcorrespondstoa pvalue shownasdottedcurvesinFigures4,5,and6fordifferentconfigu-
rationsoflineparameters.
20 Togiveaspecificexample,forthesimulationandfittingmethodsdescribed For the RS method each line spectrum was smoothed indi-
inxx3.2.1and3.2.2asetofN ¼104simulationstakes(cid:1)1dayonatop-rangePC. vidually.TheC(E)simvaluesoveranenergychannelrangeequal
Fig.6.—Broadexcess(width=0.7keV)foraspectrumcontaining800counts(left)and1600counts(right).Comparisonofthedetectionlimits,inequivalentwidth
(keV),ofthethreemethodsovertheenergybandpassofSwift.Thedataareasfollows.Dottedlines:Bayesfactoranalysis;solidlines:RSmethod;dashedlines:ppp.
Valuesbelow0.7keVhavebeenexcludedduetothewidthof thefeaturesbeinganalyzed.[SeetheelectroniceditionoftheJournalforacolorversionofthisfigure.]