Astronomy&Astrophysicsmanuscriptno.newbin cESO2016 (cid:13) January21,2016 Optimal binning of X-ray spectra and response matrix design J.S.Kaastra1,2,3 andJ.A.M.Bleeker1,3 1 SRONNetherlandsInstituteforSpaceResearch,Sorbonnelaan2,3584CAUtrecht,theNetherlands 2 LeidenObservatory,LeidenUniversity,POBox9513,2300RALeiden,theNetherlands 3 DepartmentofPhysicsandAstronomy,UniversiteitUtrecht,P.O.Box80000,3508TAUtrecht,theNetherlands January21,2016 6 ABSTRACT 1 0 Context. 2 Aims.AtheoreticalframeworkisdevelopedtoestimatetheoptimalbinningofX-rayspectra. n Methods.Wederived expressions for theoptimal bin sizefor model spectra aswellas for observed data using different levelsof a sophistication. J Results.Itisshownthatbytakingintoaccountboththenumberofphotonsinagivenspectralmodelbinandtheiraverageenergy 0 overthebinsize,thenumberofmodelenergybinsandthesizeoftheresponsematrixcanbereducedbyafactorof10 100.The − 2 responsematrixshouldthencontaintheresponseatthebincentreaswellasitsderivativewithrespecttotheincomingphotonenergy. Weprovide practical guidelines for how toconstruct optimal energy gridsaswell ashow tostructuretheresponse matrix. Afew ] examplesarepresentedtoillustratethepresentmethods. M Conclusions. I Keywords. Instrumentation:spectrographs–Methods:dataanalysis–X-rays:general . h p -1. Introduction complexthanthoseappliedtodatafrompreviousmissions,com- o putationalefficiencyisimportanttotakeintoconsideration. r tUntil two decades ago X-ray spectra of cosmic X-ray sources Forthesereasonswediscussheretheoptimalbinningofboth swere obtained using instruments such as proportional coun- a model spectra and data. In this paper, we critically re-evaluate ters and gas scintillation proportional counters with moderate [ the conceptof response matricesand the way spectra are anal- spectralresolution(typically5–20%).Withtheintroductionof ysed.Infact,weconcludethatitisnecessarytodroptheclassical 1charge-coupleddevices (CCDs) (ASCA, launch 1993) a major conceptofamatrix,andtouseamodifiedapproachinstead. vleapinenergyresolutionhasbeenachieved(upto2%)andvery 9 The outline of this paper is as follows. We start with a highresolutionhasbecomeavailablethroughgratingspectrome- 0 morein-depthdiscussionregardingthemotivationforthiswork ters,firstontheExtremeUltravioletExplorer(EUVE),andlater 3 (Sect.2).Inthefollowingsections,wediscusstheclassicalap- ontheChandraandXMM-Newtonobservatories. 5 proach to spectral modelling and its limitations (Sect. 3), the 0 These X-ray spectra are usually analysed with a forwards optimal bin size for model spectra (Sect. 4) and data (Sect. 5) 1.foldingtechnique.Firstaspectralmodelappropriatefortheob- followedbyapracticalexample(Sect.6).Wethenturntothere- served source is chosen. This model is convolved with the in- 0 sponsematrix(Sect.7),its practicalconstruction(Sect. 8), and strument response, which is represented usually by a response 6 brieflytotheproposedfileformat(Sect.9)beforereachingour matrix. The convolved spectrum is compared to the observed 1 conclusions. :spectrumandtheparametersofthemodelarevariedinafitting vprocedureinordertoobtainthebestsolution. i X This classical way of analysing X-ray spectra has been 2. Motivationforthiswork widelyadoptedandisimplemented,e.g.inspectralfittingpack- r aages such as XSPEC (Arnaud 1996), SHERPA (Freemanetal. Inthis sectionwe presentin moredepththe argumentsleading totheproposedbinningofmodelspectraandobservationaldata 2001), and SPEX (Kaastraetal. 1996). However, the applica- tion of standard concepts, such as a response matrix, is not at and our choice for the response matrix design. We do this by alltrivialforhigh-resolutioninstruments.Forexample,withthe addressingthefollowingquestions. RGS of XMM-Newton,the properlybinnedresponsematrix is 120Megabytesinsize,countingonlynon-zeroelements.Taking 2.1.Whynotusestraightforwarddeconvolution? into account that usually data from both RGS detectors and of twospectralordersarefitsimultaneously,makesitatbestslow Inhigh-resolutionopticalspectratheinstrumentalbroadeningis tohandleevenbymostpresentdaycomputersystems.Alsothe oftensmallcomparedtotheintrinsiclinewidths.Inthosecases higher spectral resolution considerably enhances the computa- itiscommonpracticetoobtainthesourcespectrumbydividing tiontimeneededtoevaluatethespectralmodels.Sincethemod- the observed spectrum at each energy by the nominal effective els appliedto Chandraand XMM-Newtondata are muchmore area(straightforwarddeconvolution). Althoughstraightforwarddeconvolution,dueto itssimplic- Sendoffprintrequeststo:J.S.Kaastra ity would seem to be attractive for high-resolutionX-ray spec- Articlenumber,page1of16 A&Aproofs:manuscriptno.newbin troscopy,it fails in severalsituations. For example,spectralor- otherandcorrespondingtothesametransition,fromnarrow-line dersmayoverlapaswiththeEUVEspectrometers(Welshetal. region lines (width few 100 kms 1), intermediate width lines − 1990) or the Chandra (Weisskopfetal. 1996) Low-Energy (1000kms 1),broad-lineregionlines(several1000uptotensof − Transmission Grating Spectrometer (LETGS; Brinkmanetal. thousandskms 1)uptorelativisticallybroadenedlines.Insuch − 2000)whentheHRC-Sdetectorisused.Inthesecasesonlycare- casessplittinginnarrow-andbroadbandfeaturesisimpossible, fulinstrumentcalibrationincombinationwithpropermodelling asthereisspectralstructureatmanydifferent(Doppler)scales. of the short wavelength spectrum can help. In case of the Re- As already stated above,there is a need for physicallyrelevant flectionGratingSpectrometer(RGS;DenHerderetal.2001)on models,regrettably,δ-functionsdonotqualifyforthis. board XMM-Newton (Jansenetal. 2001) 30% of all line flux In principle, the user could design an input model energy iscontainedinbroadwingsasaresultofscatteringonthemir- gridwherethereis onlysubstantiallyhigherspectralresolution rorandgratings,andthisfluxispracticallyimpossibletorecover at the energies where the model would predict narrow spectral bystraightforwarddeconvolution.Finally,high-resolutionX-ray features,forinstancenearforegroundabsorptionedgesorstrong spectraoftensufferfrombothseverelineblendingandrelatively sharp emission lines. However, the real source spectrum may largestatisticalerrorsbecauseoflownumbersofcountsinsome containadditionalnarrowfeaturesthatarenotanticipatedbythe partsofthespectrum.Thisrendersthemethodunsuitable. user,andwe wanttoavoidrecreatinggridsandmatricesmulti- ple times. For thisreason, the binningscheme proposedin this paperonlydependson the propertiesof the instrumentand the 2.2.Whyarehigh-resolutionspectramuchmorecomplexto observedspectrumintermsofcountsperresolutionelement. handlethanlow-resolutionspectra? The enhanced spectral resolution and sensitivity of the current 2.4.Whynotuseprecalculatedmodelgrids? high-resolution X-ray spectrometers require much more com- plexsourcemodelswithamultitudeoffreeparametersascom- Some spectral models are computationally intensive. We are paredtothecrudemodelsthatsufficeforlow-resolutionspectra awareofthemethodofprecalculatinggridsofmodels,suchthat tocoverallrelevantspectraldetailsthatcanbeexposedbythefar thespectralfittingproceedsfaster.Thisiscommonpracticewith superiorresolvingpowerofstate-of-the-artcosmicX-rayspec- for instance the APEC modelfor collisional plasmas as imple- trometers.Thisdoesnotmerelyinvolveaddingmorelinestothe mentedinXSPECorXSTARmodelsforphoto-ionisedplasmas. old models with the same number of parameters. Whereas in- For a limited set of parameters this is an acceptable solution. vestigatorsfirst fit single-temperature,solar abundancespectra, However,as hasbeenoutlinedabovein manysituations, astro- now multi-temperature, free abundance models are to be em- physicalmodelsrequiresubstantiallymorefreeparametersthan ployed for stars and clusters of galaxies, among others. Mod- canbeprovidedwithgridsfor2–4freeparameters. elsofAGNhaveevolvedfromsimplepowerlawswithaGaus- As an example, for stellar spectra not only temperaturebut sianlineprofileoccasionallysuperimposedtocomplexcontinua alsodensityandtheUVradiationfield(forHe-liketriplets)con- spectra,includingfeaturesduetoreflection,relativisticblurring, stitute importantparameters. However,the standard APEC im- softexcesses,multipleabsorbingphoto-ionisedoutflowcompo- plementation only has temperature and abundances as free pa- nents,low-ionisationemissioncontributionsoriginatingatlarge rameters. distances,whichyetagainleadstomuchlargernumbersoffree Photo-ionised,warmabsorbermodelsforAGNdependnon- parameterstobeaccountedfor. linearlyontheshapeoftheionisingspectralenergydistribution and on source parameters such as abundances and turbulence. In realistic descriptions, this requires at least five or more free 2.3.Whyisitnotpossibletomodelcomplexspectrabythe parametersandismuchtooexpensivetoconstructgrids.Using sumofasimplecontinuumplusdeltalines? gridsinsuchcasesalwaysimplieslessrealisticmodels. PracticallyallX-rayspectraintheUniversearetoocomplexto Inaddition,creatingthesegridsoffersefficiencygaininthe simulate with simple continuum shapes with superimposed δ- fittingprocedurebutshiftsthecomputationalburdentothegrid functionsmimickinglinefeatureswhenobservedathighresolu- creation. The developers of APEC communicated to us that a tion.Apartfromspectrallines,othernarrowfeaturesoftenfound fullgridcoveringallrelevanttemperaturesrequiresoftheorder are: for example, radiative recombination continua, absorption of a week computation time. A single XSTAR model (1 spec- edges, Compton shoulders, or dust features. Moreover spectral trum)takesonaverageabout20minutestocomplete,henceeven linesaretoocomplextobedealtwithbydeltafunctions,i.e.the modestgridsmaytakeaweekorsotocomplete.So,alsoitisex- astrophysicsderivingfromDopplerbroadeningornaturalbroad- tremelyuseful to be able to limit the numberof energybins in ening,or a combinationthereof(Voigtprofiles)remainstotally thosecases. unaccountedfor. Even if one would model a line emission spectrum by the 2.5.Whynotusebruteforcecomputingpower? sumofδ-lines,however,aphysicalmodelconnectingthelinein- tensitiesisneeded.Thelinefluxesaremutuallynotindependent We note that for individual spectra extensive energy grids and and,asaconsequence,manyweaklinesmaynotbedetectedin- large sizes of response matrices can be handled with present- dividuallybutaddedtogethertheymaygiveadetectablesignal. daycomputers,butthecomputationtimeinfoldingtheresponse Furthermore, except for the nearest stars, almost all X-ray matrixintothespectrumissubstantial.Thisholds,inparticular, sources are subject to foreground interstellar absorption, also ifthisprocesshastoberepeatedthousandsoftimesinspectral yielding narrow and sharp spectral features such as absorption fittinganderrorsearchesformodelswithmanyfreeparameters. edgesorabsorptionlines(forinstancethewell-knownOi1s–2p More importantly, in several cases users want to fit spectra transitionat23.5Å). oftime-variablesourcestakenatdifferentepochstogether,using As an example, Seyfert 1 galaxies show a range of emis- spectral models where some parameters are fixed between the sion lines with different widths, often superimposed on each observationswhileothersarenot.Thisrequiresloadingasmany Articlenumber,page2of16 J.S.KaastraandJ.A.M.Bleeker:Optimalbinning responsematricesasnumberofobservations,andwehaveseen convolvedwiththeinstrumentresponseR(E ,E)asfollows: ′ cases where such analyses simply cannot be carried out in this way. ∞ There are also other cases where investigators have to rely s(E′)= R(E′,E)f(E)dE, (1) on indirect methods such as making maps of line centroids or Z 0 equivalentwidths,ratherthanfullspectralfitting.Astrikingex- ample constitutes the Chandra CCD spectra of Cas A. With 1 wherethedatachannelE denotestheobserved(seebelow)pho- ′ arcsec spatial resolution the remnant covers of order 105 inde- ton energy and s has units of countss 1keV 1. The response − − pendentregions,whichshouldideallybefittogetherwithacom- functionR(E ,E)hasthedimensionsofan(effective)area,and ′ monspectralmodelwithposition-dependentparameters. canbegivenin,e.g.m2. The variable E is denoted as the observed photon energy, ′ butinpracticeitisoftensomeelectronicsignalinthedetector, 2.6.Someexplicitexamplesofspectralanalysesthatrequire the strength of which usually cannot take arbitrary values, but considerablecomputationtime can have onlya limited set of discrete values, forinstance as a A few typical examples may help to stipulate the relevance of resultofanaloguetodigitalconversion(ADC)intheprocessing ourcaseregardingrequiredcomputingtime. electronics.A goodexampleofthisisthe pulse-heightchannel Mernieretal.(2015)analysedaclusterofgalaxieswiththe foraCCDdetector.Alternatively,itmaybethepixelnumberon XMM-NewtonEPICdetectors.Thesedetectorsstillhaveamod- animagingdetectorifgratingsareused. estspectralresolution.Thespectraofthreedetectorsintwoob- Inalmostallcircumstancesitisnotpossibletocarryoutthe servationswithdifferentbackgroundlevelswerefitjointly.The integration in Eqn. (1) analytically because of the complexity model was a multi-temperature plasma model with about 50 of both the model spectrum and instrument response. For that freeparameters,includingabundancesanddifferentinstrumental reason,themodelspectrumisevaluatedatalimitedsetofener- and astrophysicalbackgroundcomponents.Fits includingerror gies,correspondingtothesamepredefinedsetofenergiesthatis searchestookabout5–10hoursonamodernworkstation.Addi- usedfortheinstrumentresponseR(E′,E).Thentheintegration tionaltimeisneededtoverifythatfitsdonotgetstuckinlocal ineqn.(1)isreplacedbyasummation.Wecallthislimitedsetof subminima.Theseauthorsareextendingtheworknowtoradial energiesthemodelenergygridor forshortthe modelgrid.For profiles in clusters, incorporating eight annuli per cluster, on a each bin j of this modelgrid, we define a lower and upperbin sample of40 clusters. Severalmonthsofcpu time usingmulti- boundaryE1j and E2j, a bin centre Ej = 0.5(E1j+E2j), and a pleworkstationsareneededtocompletethisanalysis. binwidth∆Ej = E2j E1j. − Kaastraetal. (2014a) studied the AGN Mrk 509. Chandra Providedthatthebinwidthofthemodelenergybinsissuf- HETGS spectra were fit jointly with archival XMM-Newton ficientlysmallcomparedtothespectralresolutionoftheinstru- spectra. As a result of the amount of information available by ment,thesummationapproximationto(1)isingeneralaccurate. virtue of the high spectral resolution, the spectrum was mod- TheresponsefunctionR(E′,E)hasthereforebeenreplacedbya elled using an absorption model that comprised the product of responsematrixRij,wherethefirstindexidenotesthedatachan- 28 components (four velocity components times seven ionisa- nel,andthesecondindex jthemodelenergybinnumber. tioncomponents).Withtheinclusionoffreeparametersforthe Wehaveexplicitly modellingofthecontinuumandthenarrowandbroademission features, the modelended up with about 60 free parameters. It S = R F , (2) i ij j tooktwoweeksofcomputingtimetoobtainthefinalmodelin- j cludingtheassessmentofuncertaintiesonallparameters. X Kaastraetal. (2014b) observed the Seyfert galaxy NGC where now S is the observed count rate in countss 1 for data i − 5548inanobscuredstate.StackedRGS,pn,NuSTAR,andIN- channel i and F is the model spectrum (photonsm 2s 1) for j − − TEGRAL data were analysed jointly. The model (table S3 of modelenergybin j. that paper) had 35 free parameters and consisted of six‘warm absorber’componentsandtwo‘obscured’components,givinga total of eight multiplied transmission factors, which had to be 3.2.Evaluationofthemodelspectrum determinediterativelyin abouteightsteps to reach fullconver- ThemodelspectrumF canbeevaluatedintwoways.First,the gence.Thisisbecausetheinner,obscurercomponentsaffectthe j modelcanbeevaluatedatthebincentreE ,essentiallytaking ionisingspectralenergydistributionfortheouter,warmabsorber j components. At each intermediate step a full grid of Cloudy F = f(E )∆E . (3) j j j models had to be calculated to update the transmittance of the obscurednuclearspectrum.Thefullcomputationalsotooksev- Thisisappropriateforsmoothcontinuummodelssuchasblack- eraldaysinadditiontomonthsoftrialtosetupthepropermodel bodyradiationandpowerlaws.Forline-likeemission,itismore andcalculationscheme. appropriatetointegratethelinefluxwithinthebinanalytically, taking 3. Classicalapproachtospectralmodellingandits E 2j limitations F = f(E)dE. (4) j 3.1.Responsematrices EZ 1j ThespectrumofanX-raysourceisgivenbyitsphotonspectrum f(E),afunctionofthecontinuousvariableE,thephotonenergy, Nowa seriousflawoccursinmostspectralanalysiscodes.The andhasunitsof,e.g.photonsm 2s 1keV 1.Toproducethepre- parameter S is evaluated in a straightforward way using (2). − − − i dictedcountspectrums(E )measuredbyaninstrument, f(E)is Hereby it is tacitly assumed that all photons in the model bin ′ Articlenumber,page3of16 A&Aproofs:manuscriptno.newbin j have exactly the energy E . This is necessary since all infor- puting time. Finally, the response matrices of instruments like j mation on the energy distribution within the bin is lost and F the XMM-Newton RGS become extremely large owing to ex- j isessentiallyonlythetotalnumberofphotonsinthebin.Ifthe tendedscatteringwingscausedbythegratings. model bin width ∆E is sufficiently small this is no problem, j It is thereforeimportantto keep the number of energybins however,thisis(mostoften)notthecase. assmallaspossiblewhilemaintainingtherequiredaccuracyfor An example of this is the standard response matrix for the properlinecentroiding.Fortunatelythereisamoresophisticated ASCA SIS detector (Tanakaetal. 1994), in which the investi- waytoevaluatethespectrumandconvolveitwiththeinstrument gators used a uniform model grid with a bin size of 10eV. At response. Basically, if more information on the distribution of a photon energy of 1keV, the spectral resolution (full width photonswithinabinistakenintoaccount(liketheiraverageen- at half maximum; FWHM) of the instrument was about 50eV, ergy),energygridswith broaderbins(and hence with substan- hencethelinecentroidofanisolatednarrowlinefeaturecontain- tiallyfewerbins)cangiveresultsasaccurateasfinegridswhere ingN countscanbedeterminedwithastatisticaluncertaintyof allphotonsareassumedtobeatthebincenter. 50/(2.35√N)eV. We assume here forsimplicity a Gaussian in- strumentresponse(FWHMisapproximately2.35σ,seeEq.6)). Thus, fora line with 400countsthe line centroidcan be deter- minedwithanaccuracyof1eV,tentimesbetterthanthebinsize 4. Whatistheoptimalbinsizeformodelspectra? ofthemodelgrid.Ifthetruelinecentroidisclosetothebound- aryoftheenergybin,thereisamismatch(shift)of5eVbetween 4.1.Definitionoftheproblem theobservedcountspectrumandthepredictedcountspectrumat aboutthe5σsignificancelevel.If therearemoreofthese lines In the example of the LETGS spectrum of Capella we showed inthespectrum,itispossiblethatasatisfactoryfit(e.g.accept- thattomaintainfullaccuracyforthestrongandnarrowemission able χ2 value) is never obtained, even in cases where the true lines,verysmallbinsizesarerequired,whichleadstomorethan source spectrum is known and the instrument is perfectly cali- amilliongridpointsforthemodelspectrum. brated. The problem becomes even more worrisome if, for ex- ample detailed line centroiding is performed to derive velocity Fortunately,thereareseveralwaystoimprovethesituation. fields. Firstly, the spectral resolution of most instruments is not con- A simple way to resolve these problems is just to increase stant, and one might adapt the binning to the local resolution. thenumberofmodelbins.Thisrobustmethodalwaysworks,but SomeX-rayinstrumentsusethisprocedure.Ithelps,buttheim- at the expense of a lot of computingtime. For CCD-resolution provementisnotverygoodforinstrumentslike LETGSwith a spectra this is perhaps not a problem, but with the increased differenceofonlyafactoroftwoinresolutionfromshorttolong spectral resolution and sensitivity of the grating spectrometers wavelengths. ofChandraandXMM-Newtonthisbecomescumbersome. Secondly, the finest binning is needed near features with For example, the LETGS spectrometer of Chandra large numbers of counts, for instance the strong spectral lines (Brinkmanetal. 2000) has a spectral resolution (FWHM) be- of Capella. One might therefore adjust the binning according tween0.040–0.076Åoverthe1–175Åband.A85ksobserva- to the properties of the spectrum: narrow bins near high-count tionofCapella(obsid.1248,Meweetal.2001;Nessetal.2001) regions and broader bins near low-count regions. As far as we producedN =14000countsintheFexviilineat15Å.Because know,noX-raymissionutilisesthisprocedure.Perhapsthemain line centroids can be determined with a statistical accuracy of reasonis that one first needsto knowthe observedcountspec- σ/√N,thismakesitnecessarytohaveamodelenergygridbin trum before being able to create the model energy grid for the responsematrix.Inmostcases,theredistributionmatrixiseither size σ/√N of about0.00014Å, correspondingto 278binsper precalculatedforallspectra,oriscalculatedontheflyforindi- FWHMresolutionelementandrequiring1.2millionmodelbins vidualspectra only to accountfor time-dependentinstrumental upto175Åforauniformwavelengthgrid.Althoughthissmall features. binsizeislessthanthe(thermal)widthoftheline,largermodel bin sizes would lead to a significant shift of the observed line A third way to reduce the number of bins is to revisit the with respect to the predicted line profile and a corresponding way modelspectra are calculated.As indicatedbefore,in most significant worsening of the goodnessof fit. With future, more standardprocedures,theintegratednumberofphotonsinabinis sensitiveinstrumentslikeAthena(Nandraetal.2013)suchcon- calculated.Intheconvolutionwiththeresponsematrixitisthen cernswillbecomeevenmorefrequent. tacitlyassumedthatallphotonsofthebinhavethesameenergy, Mostofthecomputingtimeinthermalplasmamodelsstems i.e. the energy Ej of the bin centre. However, it could well be fromtheevaluationofthecontinuum.Theradiativerecombina- thatthebincontainsonlyonespectralline;itmakesadifference tioncontinuumhastobeevaluatedforallenergies,forallrele- if the line is at the lower or upperlimit of the energybin or at vantionsandforallrelevantelectronicsubshellsofeachion.On the bin centre. The proper centroid of the line in the observed the other hand,the line power needsto be evaluatedonly once count spectrum, after convolution with the response matrix, is foreachline,regardlessofthenumberofenergybins.Therefore only reproducedif in the modelspectrum notonly the number thecomputingtimeisapproximatelyproportionaltothenumber of photons but also their average energy is accounted for. As ofbins.Therefore,afactorof1000increaseincomputingtime we show (see Fig. 1), accountingfor the averageenergyof the is implied from the used number of ASCA-SIS bins (1180) to photons within a bin allows us to have an order of magnitude the required number of bins for the LETGS Capella spectrum largerbinsizes. (1.2millionbins). Intheprocedureproposedherewecombineallthreeoptions Furthermore, because of the higher spectral resolution of to reduce the numberof model energybins: makinguse of the gratingspectrometerscomparedwithCCDdetectors,morecom- localspectralresolution,the strengthofthe spectralfeaturesin plexspectralmodelsareneededtoexplaintheobservedspectra, numberofphotonsN,andaccountingfortheaverageenergiesof withmorefreeparameters,whichalsoleadstoadditionalcom- thephotonswithinbins. Articlenumber,page4of16 J.S.KaastraandJ.A.M.Bleeker:Optimalbinning 4.2.Binningthemodelspectrum whereE isgivenby a 4.2.1. Zerothorderapproximation E 2j We have seen in the previoussubsection that the classical way E = f(E)EdE. (10) toevaluatemodelspectraistocalculatethetotalnumberofpho- a tonsin each modelenergybin, andto actas if all fluxis at the EZ 1j centreofthebin.Inpractice,thisimpliesthatthetruespectrum f (E) within the bin is replacedby its zeroth orderapproxima- 0 For the worst case zeroth order approximation f for f , tion, f (E),andiswrittenas 1,0 0 1,0 namely the case that f is a narrow line at the bin boundary, 0 theapproximation f yieldsexactresults.Thus,itconstitutesa f (E)=Nδ(E E ), (5) 1,1 1,0 j − majorimprovement.Infact,itiseasytoseethattheworstcase where N is the total number of photons in the bin, Ej the bin for f1,1 is a situation where f0 consists of two δ-lines of equal centreasdefinedintheprevioussection,andδistheDiracdelta- strength:oneateachbinboundary.Inthatcase,thewidthofthe function. resultingcountspectrumisbroaderthanσ.Againinthelimitof We assume, for simplicity, that the instrument response is small∆,itiseasytoshowthatthemaximumerrorδk,1for f1,1to Gaussian,centredatthetruephotonenergyandwithastandard beusedintheKolmogorov-Smirnovtestiswrittenas deviation σ. Many instruments have line spread functions (lsf) with Gaussian cores. For instrument responses with extended 1 ∆ δ = ( )2, (11) wings(e.g.aLorentzprofile)themodelbinningisalessimpor- k,1 8√2πe σ tant problem,since in the wings all spectraldetails are washed out, and only the lsf core is important. For a Gaussian profile, where e is the base of the naturallogarithms. Accordingly,the theFWHMofthelsfiswrittenas limitingbinsizefor f isexpressedas 1,1 FWHM = ln(256)σ 2.35σ. (6) ≃ ∆ 1 <2.4418λ(R)0.5N 0.25. (12) How canpwe measurethe errorintroducedwith approxima- FWHM − tion (5)? We compare the cumulative probability distribution functions(cdf)of the truespectrumandthe approximation(5), ItisseenimmediatelythatforlargeNtheapproximation f 1,1 which are both convolved with the instrumental lsf. The ap- requiresasignificantlysmallernumberofbinsthan f . 1,0 proach is described in detail in Sect. A, using a Kolmogorov- Smirnovtestwederivethemodelbinsize∆Eforwhichin97.5% ofallcasesthe approximation(5)leadstothe sameconclusion 4.2.3. Secondorderapproximation asatestusingtheexactdistribution f in testsof f versusany 0 0 Wecandecreasethenumberofbinsfurtherbynotonlycalculat- alternativemodel f . 2 ingtheaverageenergyofthephotonsinthebin(thefirstmoment The approximationeqn. (5) fails most seriously in the case ofthephotondistribution),butalsoitsvariance(thesecondmo- that the true spectrum within the bin is also a δ-function, but ment).Inthiscaseweapproximate locatedatthebinboundary,atadistance∆/2fromtheassumed linepositionatthebincentre. Themaximumdeviationδ oftheabsolutedifferenceofboth f1,2(E)= Nexp[(E Ea)/2τ2)], (13) k − cumulativedistributionfunctions,δ = F (x) F (x) (seealso k 0 1,0 eqn. A.3) occurs where f (x) = f (x|), as o−utlined a|t the end whereτisgivenby 0 1,0 of Sect. A.3. Because f (x) = f (x ∆/2), we find that the 1,0 0 − maximumoccursatx=∆/4.Aftersomealgebrawefindthatin E 2j thiscase τ2 = f(E)(E E )2dE. (14) a ∆ − δ =δ = P(∆/4) P( ∆/4)=2P(∆/4) 1= , (7) EZ k k,0 1j − − − 2√2πσ wherePisthecumulativenormaldistribution.Thisapproxima- The resulting countspectrum is then simply Gaussian with tion holdsin the limit of ∆ σ. Inserting(6) we find that the theaveragevaluecentredatE andthewidthslightlylargerthan a ≪ binsizeshouldbesmallerthan theinstrumentalwidthσ,namely √σ2+τ2. ∆ The worst case againoccurs for two δ-linesat the opposite <2.1289λkN−0.5, (8) binboundaries,butnowwithunequalstrength.Itcanbeshown FWHM inthesmallbinwidthlimitthat wherethenumberλ isdefinedby(A.9).MonteCarloresultsare k presentedinSect.4.3. 1 ∆ δ = ( )3 (15) 1,2 36√6π σ 4.2.2. Firstorderapproximation andthatthismaximumoccursforalineratioof3+√3:3 √3. A further refinement can be reached as follows. Instead of − Thelimitingbinsizefor f iswrittenas putting all photons at the bin centre, we can put them at their 1,2 averageenergy.Thisfirst-orderapproximationcanbewrittenas ∆ 2 <2.2875λ(R)1/3N 1/6. (16) − f (E)=Nδ(E E ), (9) FWHM 1,1 a − Articlenumber,page5of16 A&Aproofs:manuscriptno.newbin 4.3.MonteCarloresults Inpractice, becauseof statistical fluctuationsthe maximumfor D canalsobereachedforothervaluesofxclosetox ,andthis In addition to the analytical approximations described above, ′ m causes a somewhat more relaxed constraint on the bin size for we have performed Monte Carlo calculations as outlined in theMonteCarloresults. Sect. A.4. For our approximations of order 0, 1, and 2, corre- sponding to accounting for the number of photons only, both thenumberofphotonsandcentroidandthenumberofphotons, 4.4.Whichapproximationdowechoose? centroid, and dispersion, respectively, we found the following We nowcomparethe differentapproximations f , f , and f as approximations: 0 1 2 derivedintheprevioussubsection.Fig.1showsthattheapprox- ∆ imation f yields an order of magnitude or more improvement =min(1,y), (17) 1 FWHM overtheclassicalapproximation f0.However,theapproximation f is only slightly better than f . Moreover, the computational 2 1 with burden of approximation f is large. The evaluation of (14) is 2 0.5707 1.0 ratherstraightforward,althoughcareshouldbetakenwithsingle order0 : y= 1+ , (18) x1/2 x machineprecision;firstthe averageenergy Ea shouldbe deter- 1.404 (cid:16) 18 (cid:17) minedandthenthisvalueshouldbeusedintheestimationofτ. order1 : y= 1+ , (19) Amoreseriousproblemisthatthewidthofthelsfshouldbead- x1/4 x justedfromσto √σ2+τ2.IfthelsfisapureGaussian,thiscan 1.569(cid:16) 1.1(cid:17)4 order2 : y= 1+ , (20) becarriedoutanalytically;however,foraslightlynon-Gaussian x1/6 x1/3 lsfthetruelsfshouldbeconvolvedingeneralnumericallywith (cid:16) (cid:17) where aGaussianofwidthτtoobtaintheeffectivelsffortheparticular bin, and the computationalburden is very heavy. On the other order0 : x= N (1+0.3lnR), (21) r hand,for f onlyashiftinthelsfissufficient. 1 order1 : x= Nr(1+0.1lnR), (22) Thereforewerecommendusingthelinearapproximation f1. order2 : x= N (1+0.6lnR). (23) Theoptimalbinsizeisthusgivenby(17),(19),and(22). r Here N is the number of photons within a resolution element r andRthenumberofresolutionelements.Thepreciseuppercut- 4.5.Theeffectivearea off value of 1.0 is a little arbitrary, but we adopted it here for Above we showed how the optimal model energy grid can be simplicityasthesameasforourapproximationofthedatagrid determined,takingintoaccountthepossiblepresenceofnarrow binning(seeSect.5). spectralfeatures,thenumberofresolutionelements,andtheflux ofthesource.Wealsoneedtoaccountfortheenergydependence of the effective area, however.In the previoussection, we con- 1 sideredmerelytheeffectofthespectralredistribution(rmf);here weconsidertheeffectivearea(arf). 0.1 If the effective area Aj(E) would be a constant Aj within a modelbin j,thenforaphotonflux F inthebinthetotalcount j M 01 rateproducedbythisbinwouldbesimplyAjFj.Thisapproach H 0. isactuallyusedintheclassicalwayofanalysingspectra.Ingen- W F eral A (E) is not constant, however,and the above approach is ∆ / 0−3 justifiejdonlywhenallphotonsofthemodelbinhavetheenergy 1 of the bin centre. It is better to take into account not only the fluxF butalsotheaverageenergyE ofthephotonswithinthe 0−4 modeljbin.ThisaverageenergyE isaingeneralnotequaltothe 1 a bin centre E , and hence we need to evaluate the effectivearea j 0−5 notatEjbutatEa. 11 10 100 1000 104 105 106 107 108 109 1010 We consider here the most natural first-order extension, Nr namelytheassumptionthattheeffectiveareawithinamodelbin Fig. 1. Optimal bin size ∆ for model spectrum binning with a Gaus- isa linearfunctionoftheenergy.Foreachmodelbin j, we de- sianlsf,as afunctionof the number of counts per resolution element veloptheeffectiveareaA(E)asaTaylorseriesinE E ,which j Nr.Dash-dottedline:zerothorderapproximation(18);thicksolidline: iswrittenas − firstorderapproximation(19);dashedline:secondorderapproximation (20).Forcomparisonwealsoshowasthedottedcurvetheapproxima- A (E)=A(E )+A(E )(E E )+.... (24) j j ′ j j tion(36)forthedatagridbinning.AllresultsshownareforR=1. − The maximum relative deviation ǫ from this approximation max WeshowtheseapproximationsinFig.1.Theasymptoticbe- occurswhenEisatoneofthebinboundaries.Itisgivenby haviour is as described by (8), (12), and (16), i.e. proportional 1 toNr−1/2,Nr−1/4,andNr−1/6,respectively,althoughthenormalisa- ǫmax = 8(∆E)2A′′(Ej)/A(Ej), (25) tions for the Monte Carlo results as comparedto the analytical approximationsarehigherbyfactorsof2.27,1.63and1.38,re- where∆Eisthemodelbinwidth.Thereforebyusingtheapprox- spectively.Thedifferenceiscausedbythefactthattheanalytical imation(24)wemakeatmostarelativeerrorintheeffectivearea approximationassumedthatthemaximumoftheKolmogorov- givenby(25).Thiscanbetranslateddirectlyintoanerrorinthe Smirnovstatistic D is reachedatthe x-value x where the cu- predictedcountrate bymultiplyingǫ bythe photonflux F . ′ m max j mulativedistributionsF andF reachtheirmaximumdistance. The relative error in the count rate is thus also given by ǫ , 0 1 max Articlenumber,page6of16 J.S.KaastraandJ.A.M.Bleeker:Optimalbinning which should be sufficiently small compared to the Poissonian wherew andw arethevaluesof∆/FWHMascalculatedusing 1 a fluctuations1/√N intherelevantrange. (19)and(30),respectively. r We can do this in a more formalway by finding for which Thischoiceofmodelbinningensuresthatnosignificanter- valueofǫ atestofthehypothesisthatN isdrawnfromaPois- rorsaremadeeitherduetoinaccuraciesinthemodelortheef- r sonian distribution with average value µ is effectively indistin- fectiveareaforfluxdistributionswithinthemodelbinsthathave guishable from tests that N is drawn from a distribution with E , E . r a j mean µ(1+ǫ). We have outlined such a procedure in Sect. A, andinthelimitoflargeµwewrite 5. Databinning µ+q √µ=µ(1+ǫ)+q µ(1+ǫ), (26) α kα 5.1.Introduction whereαisthesizeofthetepstandqα isgivenbyG(qα) = 1−α Most X-ray detectors count the individual photons and do not with G the cumulative normal probability distribution. Further registertheexactenergyvaluebutadigitisedversionofit.Then kαisthesizeofthetestwhenweusetheapproximation(24).In a histogram is produced containing the number of events as a thelimitofsmallǫ,thesolutionof(26)isgivenby functionof theenergy.The binsize of thesedata channelside- ally should not exceed the resolution of the instrument, other- ǫ p(α,k)/ (N ), (27) r ≃ wise important information may be lost. On the other hand, if where we hapve approximated µ by the observed value Nr, and the bin size is too small, one may have to dealwith low statis- where p(α,k) = qα qkα. Forα = 0.025andk = 2, we obtain ticsperdatachannel,insufficientsensitivityofstatisticaltestsor − p(α,k)=0.31511.Combiningtheseresultswith(25),weobtain a large computationaloverhead caused by the large number of anexpressionforthemodelbinsize data channels.Low numbersofcountsper data channelcan be alleviatedbyusingC-statistics(oftenslightlymodifiedfromthe ∆ = 8Ap(α,k) 0.5 Ej N 0.25. (28) originaldefinition by Cash 1979); computationaloverheadcan FWHM E2jA′′ FWHM r− be a burden for complex models, but insufficient sensitivity of (cid:16) (cid:17) (cid:16) (cid:17) statisticaltestscanleadtoinefficientuseoftheinformationthat Thebin widthconstraintderivedhere dependsuponthe di- is contained in a spectrum. We illustrate this inefficient use of mensionless curvature of the effective area A/E2A . In most j ′′ informationinAppendixC.Inthissectionwederivetheoptimal partsoftheenergyrangethisisanumberoforderunityorless. binsizeforthedatachannels. SincethesecondprefactorE /FWHMisbydefinitiontheresolu- j tionoftheinstrument,weseebycomparing(28)with(12)that, in general, (12) gives the most severe constraint upon the bin 5.2.TheShannontheorem width. This is the case unless either the resolution becomes of TheShannon(1949)samplingtheoremstatesthefollowing:Let orderunity,whichhappens,forexamplefortheRosat(Truemper f(x) be a continuoussignal. Let g(ω) be its Fourier transform, 1982) PSPC (Pfeffermannetal. 1987) detectorat low energies, givenby or the effective area curvature becomes large, which may hap- pen,forexampleneartheexponentialcut-offscausedbyfilters. ∞ Large effective area curvature due to the presence of expo- iωx g(ω)= f(x)e dx. (32) nentialcut-offsisusuallynotaproblem,sincethesecut-offsalso Z causethecountratetobelowandhenceweakenthebinningre- −∞ quirements.Ofcourse,discreteedgesintheeffectiveareashould Ifg(ω)=0forall ω >W foragivenfrequencyW,then f(x)is alwaysbeavoidedinthesensethatedgesshouldalwayscoincide bandlimited,andi|nt|hatcaseShannonhasshownthat withbinboundaries. Inpractice,itisalittlecomplicatedtoestimatefrom,forex- ∞ sinπ(x/∆ n) amplealook-uptableoftheeffectiveareaitscurvature,although f(x)= fs(x)≡ f(n∆) π(x/∆ −n) . (33) thisisnotimpossible.Asasimplificationfororderofmagnitude nX=−∞ − estimates, we use the case where A(E) = A0ebE locally, which In (33), the bin size ∆ = 1/2W. Thus, a band-limitedsignal is afterdifferentiationyields completely determined by its values at an equally spaced grid withspacing∆. 8A dlnE = √8 . (29) The above is used for continuoussignals sampled at a dis- E2A dlnA cretesetofintervalsx.However,X-rayspectraareessentiallya s j ′′ i histogramofthenumberofeventsasafunctionofchannelnum- Insertingthisinto(28),weobtainourrecommendedbinsize,as ber.Wedonotmeasurethesignalatthedatachannelboundaries, faraseffectiveareacurvatureisconcerned butwemeasurethesum(integral)ofthesignalbetweenthedata ∆ dlnE E channelboundaries.HenceforX-rayspectraitismoreappropri- FWHM =1.5877 dlnA FWHjM Nr−0.25. (30) atetostudytheintegralof f(x)insteadof f(x)itself. Wescale f(x)torepresentatrueprobabilitydistribution.The (cid:16) (cid:17)(cid:16) (cid:17) cumulativeprobabilitydensitydistributionfunctionF(x)iswrit- 4.6.Finalremarks tenas In the previous two subsections we have given the constraints x fordeterminingtheoptimalmodelenergygrid.Combiningboth F(x)= f(y)dy. (34) requirements(19)and(30)weobtainthefollowingoptimalbin Z size: −∞ ∆ 1 IfweinserttheShannonreconstruction(33)in(34),afterinter- = , (31) FWHM 1/w +1/w changingtheintegrationandsummationandkeepingintomind 1 a Articlenumber,page7of16 A&Aproofs:manuscriptno.newbin that we cannot evaluate F(x) at all arbitrary points but only at plottedinFig.3.WehavealsoverifiedthiswiththeMonteCarlo thosegridpointsm∆forintegermwherealso f issampled,we methoddescribedbefore. s obtain ∆ ∞ π F (m∆)= f(n∆) +Si[π(m n)] . (35) s π 2 − 1 n= X−∞ n o ThefunctionSi(x)isthesine-integralasdefined,forexamplein 0.8 Abramowitz&Stegun(1965).Theexpression(35)forF equals s M F if f(x) is band limited. In that case at the grid points F is H W completelydeterminedby the value of f at the grid points. By F invertingthisrelation,onecouldexpress f atthegridpointsas ∆/0.6 a unique linear combination of the F-values at the grid. Since Shannon’stheoremstatesthat f(x)forarbitraryxisdetermined 0.4 completelybythe f-valuesatthegrid,weinferthat f(x)canbe completelyreconstructedfromthediscretesetofF-values.And then, by integrating this reconstructed f(x), F(x) is also deter- 0.2 mined. 1 10 100 1000 104 105 106 107 108 109 1010 1011 1012 WeconcludethatF(x)isalsocompletelydeterminedbythe Nr setofdiscretevaluesF(m∆)at x = m∆forintegervaluesofm, Fig.3.Optimalbinsize∆fordatabinningwithaGaussianlsfasafunc- providedthat f(x)isbandlimited. tionofthenumber ofcountsperresolutionelement N.Dottedcurve: r Fornon-band-limitedresponses,weuse(35)toapproximate analyticalapproximationusingSect.A.3;stars:resultsofMonteCarlo the true cumulative distribution function at the energy grid. In calculation;solidline:ourfinallyadoptedbinsizeEq.(36)basedona doing this, a small error is made. The errors can be calculated fittoourMonteCarloresultsforlargevaluesofNr;dashedline:com- monlyadoptedbinsize1/3FWHM. easily by comparing F (m∆) with the true F(m∆) values. The s binning∆issufficientifthesamplingerrorsaresufficientlysmall compared with the Poissonian noise. We elaborate on what is ItisseenthattheMonteCarloresultsindeedconvergetothe sufficientbelow. analytical solution for large Nr. For Nr = 10, the difference is 20%.WehavemadeasimpleapproximationtoourMonteCarlo results(shownasthesolidlineinFig.3),wherewecutitofftoa 5.3.Optimalbinningofdata constantvalueforN oflessthanabout2.BoththeMonteCarlo r resultsandtheanalyticalapproximationbreakdownatthissmall numberofcounts.Infact,forasmallnumberofcountsbinning 1 toabout1FWHMissufficient. 0.1 We alsodeterminedthedependenceonthenumberofreso- 0.01 lutionelementsRforvaluesofRof1,10,100,1000,and10000. 10−3 Ingeneral,theresultsevenforR=104arenottoodifferentfrom 10−4 the R = 1 case. We founda verysimple approximationfor the ∆m)|10−5 dependenceon R that is well described by the following func- F(10−6 tion: − ∆)10−7 m max |F(s11110000−−−−119810 FW∆HM = 10.08+17+.05/.x9+/x1.8/x2 iiffxx≤>22..111199;, (36) 10−12 with  10−13 10−14 0.2 0.5 1 2 x ln[Nr(1+0.20lnR)]. (37) ∆/FWHM ≡ Fig.2.Maximumdifferenceforthecumulativedistributionfunctionof We note that for digitisation of analogue electronic signals aGaussiananditsShannonapproximation,asafunctionofthebinsize that represent spectral information, a binning (=channel) cri- ∆.Themaximumisalsotakenoverallphasesoftheenergygridwith terion of 1/3 FWHM or better is often adopted (e.g. Davelaar respecttothecentreoftheGaussian. 1979),sinceonedoesnotknowapriorithelevelofsignificance ofthemeasuredsignal(seethedashedlineinFig.3). WeapplythetheoryoutlinedabovetothecaseofaGaussian redistributionfunction.We firstdeterminethemaximumdiffer- 5.4.Finalremarks ence δ between the cumulative Gaussian distribution function k F0(x)anditsShannonapproximationF1(x) = Fs(x),asdefined We have estimated conservative upper bounds for the required by (A.3). We show this quantity in Fig. 2 as a function of the databinsize.Inthecaseofmultipleresolutionelements,wehave binsize∆.ItisseenthattheShannonapproximationconverges determinedtheboundsfortheworstcasephaseofthegridwith quickly to the true distribution for decreasing values ∆. Using respecttothedata.Inpractice,itisnotlikelythatallresolution δ =λ /√N(see(A.9)),foranygivenvalueofNandλ =0.122 elementswouldhavetheworstpossiblealignment.However,we k k k wecaninverttherelationandfindtheoptimalbinsize∆.Thisis recommendusingtheconservativeupperlimitsasgivenby(37). Articlenumber,page8of16 J.S.KaastraandJ.A.M.Bleeker:Optimalbinning Another issue is the determination of N , the number of 7. Theresponsematrix r eventswithintheresolutionelement.Wehavearguedthatanup- Wehaveshownintheprevioussectionshowtheoptimalmodel perlimittoN canbeobtainedbycountingthenumberofevents r anddataenergygridscanbeconstructed.Wehaveproposedthat within one FWHM and multiplying it by the ratio of the total for the evaluation of the model spectrum both the number of areaunderthelsf(shouldbeequalto1)totheareaunderthelsf photons in each bin as well as their average energy should be withinoneFWHM(shouldbesmallerthan1).FortheGaussian determined.Wenowdeterminehowthisimpactstheconceptof lsf,thisratioequals1.314;forotherlsfs,itisbettertodetermine theresponsematrix. therationumericallyandnottousetheGaussianvalue1.314. In order to acquire high accuracy,we need to convolvethe For the Gaussian lsf, the resolution depends only weakly modelspectrum for the bin, approximatedas a δ-functioncen- upon the number of counts N (see Fig. 3). For low count rate r tredaroundE ,withtheinstrumentresponse.Inmostcaseswe partsofthespectrum,thebinningruleof1/3FWHMusuallyis a cannotdothisconvolutionanalytically,sowehavetomakeap- tooconservative. proximations.Fromourexpressionsfortheobservedcountspec- Finally, in some cases the lsf may contain more than one trum s(E ), eqns. (1) and (2), it can be easily derived that the component.For example,the grating spectra have higherorder ′ numberofcountsorcountratefordatachanneliisgivenby contributions. Other instruments have escape peaks or fluores- centcontributions.Ingeneralitisadvisabletodeterminethebin E size foreachofthesecomponentsindividuallyandsimplytake i′2 ∞ thesmallestbinsizeasthefinalone. S = dE R(E ,E)f(E)dE, (38) i ′ ′ Z Z E 0 6. Apracticalexample i′1 where, as before, E and E are the formal channel limits for In order to demonstrate the benefits of the proposed binning i′1 i′2 data channeli and S is the observedcountrate in counts/sfor schemes,Wehaveappliedthemtothe85ksChandraLETGSob- i data channel i. Interchanging the order of the integrations and servationofCapella mentionedearlier. Thespectrumspansthe definingthemono-energeticresponsefordatachannelibyR˜ (E) 0.86–175.55Å rangewith a spectralresolution(FWHM) rang- i asfollows: ingbetween0.040–0.076Å.InTable1weshowthenumberof databins,modelbins,andresponsematrixelementsfordifferent E i′2 rebinningschemesofthisspectrum.Forsimplicityweconsider R˜ (E) R(E ,E)dE , (39) onlythepositiveandnegativefirst-orderspectrum,andthatthe i ′ ′ ≡ non-zero elements of the response matrix span a full range of Z E fourtimestheinstrumentalFWHM.Weapplythealgorithmsas i′1 described in this paper to generate the data and model energy weobtain grids. The highest resolution of 0.040 Å occurs at short wave- ∞ lengths, and the strongest line is the Fexvii blendat 17.06and Si = f(E)R˜i(E)dE. (40) 17.10Å,withamaximumnumberofcountsN of15000counts. Z r 0 The first case we consider (case A afterwards) is that we adopta data andmodelgridwith constantstep size forthe full From the above equation we see that as long as we are inter- range. Accounting for the maximum Nr value, we have a data estedintheobservedcountrateSi ofagivendatachanneli,we get that number by integrating the model spectrum multiplied binsizeof0.02Åandamodelbinsizeof0.14mÅ. by the effective area R˜ (E) for that particular data channel. We Case B is similar to case A except that we account for the i haveapproximated f(E)foreachmodelbin jby(9),sothat(40) variable resolution of the instrument. This only decreases the becomes numberofbinsbyafactorof1.375. Incase Cwe accountforthenumberofcountsin eachres- olution element, but we still keep the ”classical” approach of Si = FjR˜i(Ea,j), (41) puttingallphotonsatthecentreofthemodelbins.Thenumber j X ofresolutionelementsRinthespectrumis3110. whereasbeforeE istheaverageenergyofthephotonsinbin IncaseDweaccountfortheaverageenergyofthephotons a,j jgivenby(10),and F isthetotalphotonfluxforbin j,ine.g. withinabin.Thenumberofmodelbinsandresponseelements j neededdropsbymorethananorderofmagnitudeforthiscase, photonsm−2s−1.Itisseenfrom(41)thatweneedtoevaluateR˜i notatthebincentreE butatE ,asexpected. ascomparedtocaseC. j a,j FormallywemaysplitupR˜ (E)inaneffectiveareapartA(E) i Table 1. Number of data bins, model bins, and response matrix ele- andaredistributionpartr˜(E)insuchawaythat i mentsfordifferentrebinningschemesoftheChandraLETGSspectrum ofCapella. R˜ (E)= A(E)r˜(E). (42) i i Binning Databins Modelbins Response Wehavechosenourbinningalreadyinsuchawaythatwehave elements sufficientaccuracywhenthetotaleffectiveareaA(E)withineach A 8.52 103 1.26 106 1.01 107 modelenergygridbin jisapproximatedbyalinearfunctionof B 6.20×103 9.14×105 7.31×106 thephotonenergy E.Hencethearf-partofR˜i isofnoconcern. × × × We only need to check how the redistribution(rmf) part r˜ can C 5.12 103 1.23 105 9.84 105 i × × × becalculatedwithsufficientlyaccuracy. D 5.12 103 8.21 103 6.57 104 × × × Forr˜i theargumentsareexactlythesameasforA(E)inthe sensethatifweapproximateitlocallyforeachbin jbyalinear Articlenumber,page9of16 A&Aproofs:manuscriptno.newbin functionofenergy,themaximumerrorthatwemakeispropor- 8. Constructingtheresponsematrix tionaltothesecondderivativewithrespecttoEofr˜,cf.(25). i Intheprevioussectionweoutlinedthebasicstructureofthere- Infact, for a Gaussian redistributionfunctionthe following sponsematrix. Foran accuratedescriptionof the spectrum,we isstraightforwardtoprove: needboththeresponsematrixanditsderivativewithrespectto Theorem1.Assumethatforagivenmodelenergybin jallpho- thephotonenergy.Webuildonthistoconstructgeneralresponse tons are located at the upperbin boundary E +∆/2. Suppose matricesformorecomplexsituations. j thatforalldatachannelsweapproximater˜ byalinearfunction i of E, and the coefficients are the first two terms in the Taylor 8.1.Dividingtheresponseintocomponents expansionaroundthebincentre E .Thenthemaximumerrorδ j made in the cumulative countdistribution (as a function of the Usuallyonlythenon-zeromatrixelementsofaresponsematrix datachannel)isgivenby(12)inthelimitofsmall∆. arestoredandused.Thisisdonetosavebothstoragespaceand computationaltime.Theprocedureas usedinXSPEC (Arnaud Theimportanceoftheabovetheoremisthatitshowsthatthe 1996) and the older versions of SPEX (Kaastraetal. 1996) is binningforthemodelenergygridthatwehavechoseninSect.4 thatforeach modelenergybin j the relevantcolumnof there- is also sufficiently accurate so that r˜i(E) can be approximated sponse matrix is subdividedinto groups.A groupis a contigu- by a linear functionof energywithina modelenergybin j, for ouspiece of the column with non-zeroresponse. These groups eachdatachanneli.Sincewealreadyshowedthatourbinningis are stored in a specific order; starting from the lowest energy, alsosufficientforasimilarlinearapproximationtoA(E),italso allgroupsbelongingtoasinglephotonenergyarestoredbefore followsthat the totalresponse R˜i(E) can be approximatedby a turningtothenexthigherphotonenergy. linearfunction.Hence,withinthebin jweuse Thisisnotoptimalneitherintermsofstoragespacenorcom- putational efficiency, as illustrated by the following example. dR˜ R˜ (E )=R˜ (E )+ i(E ) (E E ). (43) For the XMM-Newton/RGS, the response consists of a narrow i a,j i j j a,j j dEj − Gaussian-likecorewithabroadscatteringcomponentduetothe gratingsinaddition.TheFWHMofthescatteringcomponentis Inserting the above in (41) and comparing with (2) for the 10timesbroaderthan thecoreof theresponse.Asa result, if classicalresponsematrix,wefinallyobtain ∼ theresponseissavedasa classical matrix,we endupwith one responsegroupperenergy,namelythecombinedcoreandwings S = R F +R (E E )F , (44) responsebecausetheyoverlapintheregionoftheGaussian-like i ij j ′ij a,j− j j core. As a result, the response file becomes large. This is not j X necessarybecausethe scattering contributionwith its ten times whereRijistheclassicalresponsematrix,evaluatedforphotons largerwidth needsto be specified only ona modelenergygrid atthebincentre,andR′ijisitsderivativewithrespecttothepho- with ten times fewer bins, comparedto the Gaussian-like core. ton energyE . Inadditionto the classical convolution,we thus Thus,byseparatingoutthecoreandthescatteringcontribution, j get a second term containing the relative offset of the photons the totalsize of the responsematrix can be reducedbyabouta withrespecttothebincentre.Thisisexactlywhatweintended factorof10.Ofcourse,asaconsequenceeachcontributionneeds tohavewhenwearguedthatthenumberofbinscouldbereduced tocarryitsownmodelenergygridwithit. considerablyby just taking that offset into account.It is just at Thereforewe propose subdividingthe response matrix into theexpenseofanadditionalderivativematrix,whichmeansonly itsconstituentcomponents.Thenforeachresponsecomponent, afactoroftwomorestoragespaceandcomputationtime.Butfor the optimal model energy grid can be determined according to this extra expenditurewe gainedmuch morestorage space and themethodsdescribedinSect.4,andthismodelenergygridfor computational efficiency because the number of model bins is the componentcan be stored togetherwith the responsematrix reducedbyafactorbetween10–100. partofthatcomponent.Furthermore,atanyenergyeachcompo- Finally we make a practicalnote. The derivativeR can be nentmay have at mostone response group.If there were more ′ij calculatedinpracticeeitheranalyticallyorbynumericaldiffer- responsegroups,thecomponentshouldbesubdividedfurther. entiation. In the last case, it is more accurate to evaluate the InX-raydetectorsotherthanthe RGSdetectorthesubdivi- derivativebytakingthedifferenceatE +∆/2andE ∆/2,and, sion could be different. For example, with CCD detectors one j j − whereverpossible,notto evaluateitat oneof these boundaries could split up the response into for components: the main di- andthebincentre.Thislastsituationisperhapsonlyunavoidable agonal, the Si fluorescence line, the escape peak, and a broad atthefirstandlastenergyvalue. componentduetosplitevents. Also, negative response values should be avoided. Thus it shouldbeensuredthatR +R hisnon-negativeeverywherefor ij ′ij 8.2.Complexconfigurations ∆/2 h ∆/2.Thiscanbetranslatedintotheconstraintthat − ≤ ≤ R shouldbelimitedalwaystothefollowinginterval: In most cases an observer analyses the data of a single source ′ij with a single spectrum and response for a single instrument. 2R /∆E R 2R /∆E. (45) However,morecomplicatedsituationsmayarise.Examplesare: − ij ≤ ′ij ≤ ij Whenever the calculated value of R should exceed the above 1. A spatially extended source, such as a cluster of galaxies, ′ij limits, thelimitingvalueshouldbeinsertedinstead.Thissitua- withobservedspectraextractedfromdifferentregionsofthe tion may happen,for example for a Gaussian redistributionfor detector,butwiththeneedtobeanalysedsimultaneouslydue responsesafewσawayfromthecentre,wheretheresponsefalls totheoverlapinpoint-spreadfunctionfromoneregiontothe off exponentially. However, the response R is small for those other. ij energies anyway, so this limitation is not serious; this is only 2. Forthe RGSofXMM-Newton,theactualdataspacein the becausewewanttoavoidnegativepredictedcountrates. dispersiondirectionisactuallytwo-dimensional:theposition Articlenumber,page10of16