An Unfolding Method for X-ray Spectro-Polarimetry
Fabian Kislat, Matthias Beilicke, Qingzhen Guo, Anna Zajczyk, Henric Krawczynski
AAn Unfolding Method for X-ray Spectro-Polarimetry
F. Kislat ∗ , M. Beilicke, Q. Guo, A. Zajczyk, H. Krawczynski Washington University in St. Louis, Department of Physics and McDonnell Center for the Space Sciences, One Brookings Dr., CB 1105,St. Louis, MO 63130
Abstract
X-ray polarimetry has great scientific potential and new experiments, such as X-Calibur, PoGOLite, XIPE, and GEMS,will not only be orders of magnitude more sensitive than previous missions, but also provide the capability to measurepolarization over a wide energy range. However, the measured spectra depend on the collection area, detector responses,and, in case of balloon-borne experiments, the absorption of X-rays in the atmosphere, all of which are energy dependent.Combined with the typically steep source spectra, this leads to significant biases that need to be taken into account tocorrectly reconstruct energy-resolved polarization properties. In this paper, we present a method based on an iterativeunfolding algorithm that makes it possible to simultaneously reconstruct the energy spectrum and the polarizationproperties as a function of true photon energy. We apply the method to a simulated X-Calibur data set and show thatit is able to recover both the energy spectrum and the energy-dependent polarization fraction.
Keywords:
X-rays, Polarization, Unfolding, X-Calibur
1. Introduction
X-ray polarimetry holds the promise to resolve the in-ner regions of compact systems like mass accreting blackholes in X-ray binaries and X-ray bright neutron stars [1].For example, spectropolarimetric observations of pulsarsand pulsar wind nebulae can constrain the geometry andlocale of particle acceleration in these sources. Measure-ments of the polarization of X-rays from the Crab Neb-ula indicate an increase in the polarization fraction and achange in polarization angle in the energy range betweena few keV and
100 keV indicating that the γ -ray emissionmust come from a small, highly ordered region, whereasX-rays are emitted from all morphological features of thepulsar wind nebula. Spectropolarimetric observations canconstrain the magnetic structure of jets in Gamma RayBursts and Active Galactic Nuclei. Furthermore, X-raypolarimetry can be used to measure the masses and spinsof black holes and the orientation of their inner accretiondisk, as well as accretion disks and accretion disk coro-nae of Active Galactic Nuclei. See Ref. [1] and referencestherein for more details.While this potential has long been appreciated, theOSO-8 satellite launched in 1978 has been the only missionwith a dedicated X-ray polarimeter so far that measuredX-ray polarization of an astrophysical source [2]. Newtechnological developments enabled the design of compactwide-bandpass polarimeters with a large collection areasuch as the proposed satellites GEMS [3] and XIPE [4] or ∗ Corresponding author
Email address: [email protected] (F. Kislat) the balloon-borne hard X-ray polarimeters X-Calibur [5]and PoGOLite [6]. The collection areas and detection effi-ciencies of these experiments will allow spectropolarimetricobservations with unprecedented sensitivity. In this paperwe study statistical methods that can be used to analyzethe data of such experiments.The measurement of the polarization fraction and di-rection of the abovementioned experiments make use ofthe photo-electric effect (GEMS, XIPE) or the Comptoneffect (PoGOLite, X-Calibur).Photoelectrons are emitted preferentially parallel tothe electric field of the electromagnetic wave associatedwith the photon. The differential cross section of the pho-toelectric effect in the non-relativistic case can be approx-imated as [7]: d σ dΩ = r Z α (cid:18) m e c E (cid:19) √ θ cos φ (1 − β cos θ ) , (1)where r is the classical electron radius, Z the atomic num-ber of the absorbing material, α the fine-structure con-stant, m e the electron rest mass, E the photon energy, θ the angle between the incoming photon and the emittedphotoelectron, β its speed in units of c , and φ the azimuthangle of the emitted electron with respect to the polariza-tion direction of the incident X-ray.Photons scatter preferentially perpendicular to the elec-tric field, as governed by the Klein-Nishina cross section(see e. g. Ref. [8]): d σ dΩ = r k k (cid:20) k k + k k − θ cos η (cid:21) , (2) Preprint submitted to Astroparticle Physics October 3, 2018 a r X i v : . [ a s t r o - ph . I M ] N ov here η is the angle between the electric vector of theincident photon and the scattering plane, k and k arethe wave-vectors before and after scattering, and θ is thescattering angle.The azimuthal scattering angle or emission directionof the photoelectron is, therefore, a proxy for the polar-ization angle. In general, scattering polarimeters measurethe azimuthal scattering angle of an X-ray photon and theenergy of the scattered photon. In some realizations anX-ray mirror is used to focus the beam onto the detectorassembly.Since typical source spectra exhibit steep power-laws,the energy resolution and energy-dependent detection effi-ciency have to be taken into account. The most importanteffects that must be considered are:• the energy resolution of the detector;• the energy lost in the scattering process, which is notmeasured by all experiments;• the energy-dependent effective area of grazing inci-dence mirrors;• absorption of photons in the atmosphere in case ofballoon-borne instruments;• the energy-dependent detection efficiency.In this paper we describe an unfolding algorithm, andshow that it can be used to determine flux and polarizationfraction and direction as a function of photon energy withsmall biases.In Section 2 we will define the problem. In Section 3we will introduce an unfolding method that can be usedto reconstruct the photon spectrum of the source whilepreserving the energy-dependent azimuth distribution ofevents. In Section 4, we apply the method to a set of sim-ulated X-Calibur data. Finally, in Section 5, we summarizeour findings and present our conclusions.
2. Formulation of the problem
Typical X-ray polarimeters measure the energy of thescattered photon or the energy of the recoil electron. Thisenergy will differ from the energy of the incident photon byan unknown amount governed by underlying physical pro-cess (e. g. photoelectron emission or Compton scattering)and the finite energy resolution of the detectors.Grazing-incidence X-ray mirrors focus the X-rays froma source onto a detector and make it possible to combinelarge detection areas with rather small detectors, and thusto achieve excellent signal to background ratios. In addi-tion, they change the polarization properties of the inci-dent photons by < [9, 10]. However, their effective ar-eas depend strongly on the energies of incident photons. Incase of a balloon-borne experiment, the atmosphere abovethe detector absorbs low-energy photons. For illustration,Fig. 1 shows the effective collection area of the InFOC µ S Photon energy / keV9 10 20 30 40 50 60 70 E ff e c t i v e a r ea / c m ° Zenith 0 ° Zenith 30 ° Zenith 60 → mirror dominated atmosphere dominated ← Figure 1: Effective collection area of a balloon-borne polarimeterflown at an altitude of ∼
45 km (atmospheric depth of . − at ◦ elevation), taking into account the absorption of X-rays in theresidual atmosphere. At low energies the limiting factor is absorptionin the atmosphere, while at high energies the mirror effective arealimits the collection area. mirror [11], folded with the energy-dependent transmis-sivity of the atmosphere (from [12]) for a flight altitudeof
45 km ( . − at ◦ elevation).Together, the strong energy dependence of the effectivecollection area, the finite energy resolution, and the steeplyfalling source spectra, lead to a significant distortion of themeasured energy spectra, which needs to be taken intoaccount in a spectropolarimetric analysis.Mathematically, these effects can be described by aconvolution of a true spectrum with a detector responsefunction which includes effects of the mirror and the at-mosphere. If the spectrum is measured in n e discrete en-ergy bins, the measurement can be expressed as a vector ofevent counts N ei ( i = 1 , . . . , n e ) and the detector responseis represented by the response matrix R , which relates thetrue number of photons N cj ( j = 1 , . . . , n c ) in the j -thenergy bin to the observed counts in the i -the bin: N ei = n c (cid:88) j =1 R ij N cj . (3)Each entry in the response matrix corresponds to the prob-ability to measure an energy in bin i given a true energyin bin j : R ij = P ( E ei | E cj ) . (4)More generally, what was refered to as “true spectrum”so far is a set of “ causes ”, Φ cj , with a set of properties –one of them being the energy of the incident photon –,which lead to a set of “ effects ”, Φ ei , with a set of propertiesthat can be measured in the experimental setup (the “mea-sured spectrum”) – one of them for example the energy ofthe scattered photon. Therefore, instead of the “measuredspectrum” we will from now on refer to the measured ef- ects . The response matrix R ij then describes the proba-bility that cause j will lead to effect i .In general, neither the binning of nor the number ofproperties associated with Φ e and Φ c need to be the same.Thus, it is straight-forward to generalize equation (3) toinclude information about the azimuthal scattering angle φ in cause bins k = 1 , . . . , n cφ and effect bins (cid:96) = 1 , . . . , n eφ respectively. Then, each cause represents a combinationof true energy and scattering angle, Φ cj = ( E cp , φ ck ) , j =1 , . . . , n cE n cφ , and correspondingly each effect represents acombination of measured energy and scattering angle. Thevectors N c and N e now represent event numbers in binsof these more general, higher-dimensional causes Φ cj andeffects Φ ei .Additional observables (i. e. parameters of the effects )or parameters of the causes (i. e. output parameters of theunfolding) can be added in the same way. For instance,when analysing the X-Calibur data in Section 4, we addthe coordinate of the observed photon along the opticalaxis as an additional input parameter in order to improvethe energy resolution.The response matrix is normalized such that for a given cause bin j the sum over all effect bins is the detectionefficiency ε j for photons in the j -th cause bin: (cid:88) i R ij = ε j , (5)with ≤ ε j ≤ .In general, measurements will contain a significant frac-tion of background events caused by cosmic rays and al-bedo photons [5], despite active and passive shielding ofthe detector. If the background distributions of the inputvariables are known – ideally, they should be measuredduring flight in special off-source data taking runs –, thebackground can simply be subtracted from the input dis-tributions. This is possible because the unfolding methoddoes not consider individual events but only takes distri-butions of measured quantities as input. Therefore, in theremainder of this paper, we will neglect any backgroundand only consider signal events.According to the Klein-Nishina cross section, Eq. (2),photons scatter preferentially perpendicular to their elec-tric field vector; and according to Eq. (1) photoelectronsare emitted preferentially parallel to the polarization vec-tor. This results in a sinusoidal modulation of the az-imuthal scattering distribution with a ◦ period withminima at the direction of the polarization plane, and arelative amplitude proportional to the degree of polariza-tion. A common method to determine polarization param-eters is to fit a sine function to the azimuthal distribution.The modulation amplitude µ = C max − C min C max + C min (6)is therefore a measure of the polarization fraction. C max and C min are the maximum and minimum of the sinu- soidal modulation of the azimuthal distribution. The per-formance of a polarimeter can be characterized by themodulation obtained when observing a polarizedsource, µ , called the modulation factor .
3. Energy-resolved polarimetry
Typically, the response matrix R is poorly conditionedand simply inverting it to obtain an estimate of N c froma measured spectrum N e leads to unnatural fluctuationsof the result (see e. g. Ref. [13]). Instead, we use an itera-tive unfolding technique [14], which avoids this undesiredeffect. In the following we summarize this method in theapplication to spectropolarimetric data.Starting from a prior distribution P ( r ) (Φ cj ) in the r -thiteration, the posterior conditional probability P ( r ) (Φ cj | Φ ei ) is calculated using Bayes’ Theorem: P ( r ) (Φ cj | Φ ei ) = P (Φ ei | Φ cj ) P ( r ) (Φ cj ) (cid:80) k P (Φ ei | Φ ck ) P ( r ) (Φ ck ) . (7)Then, an estimate of the true spectrum is obtained: ˆ N c ( r ) j = 1 ε j (cid:88) i N ei P ( r ) (Φ cj | Φ ei ) , (8)where ε j is the efficiency in cause bin j as defined inEq. (5). At the end of each iteration P ( r ) (Φ cj ) is replacedby P ( r +1) (Φ cj ) = ˆ N c ( r ) j (cid:80) p ˆ N c ( r ) p . (9)This procedure is repeated for a number of iterations, n iter ,which has to be determined in advance, typically usingMonte Carlo simulation studies (see Section 3.3).Unlike suggested in Ref. [14], ˆ N c ( r ) j was not smoothedbetween iterations in the studies presented in this paper.Smoothing may result in a “nicer” looking result. How-ever, it can potentially introduce additional systematic un-certainties if the results depend on the smoothing kernelbeing used.The initial prior P (0) (Φ cj ) is chosen either based onprior knowledge of the source or on the observed spec-trum. One can either start from a flat azimuthal distri-bution or from an azimuthal distribution derived from theobserved data. Since the unfolding result might dependon the choice of prior, the influence of this choice shouldbe studied by repeating the unfolding with a range of “rea-sonable” priors.The result of this unfolding is a two-dimensional distri-bution of the best estimate of event numbers as a functionof incident photon energy and scattering angle. A pro-jection of this distribution onto the energy axis yields thesource spectrum. Slices along the azimuth axis can be fit-ted with a sine function in order to obtain polarizationfraction and angle for individual energy bins.3 .2. Determination of the response matrix In general, the response of the detector, R , can bedetermined from Monte Carlo simulations. These simula-tions have to accurately describe all parts of the experi-ment:• the absorption of X-ray photons in the atmosphere;• the effective collection area of the mirror;• the physical processes in the Compton scatteringmedium;• the response of the photon detectors;• and the resolution of the readout electronics.For most astrophysical applications, it is reasonable tosimulate a powerlaw spectrum of incident photons. Foreach simulated photon that triggers a detector, one thendetermines the observed electronic signal at the readoutelectronics. In this way, the distribution of measured ob-servables is determined for each true cause bin, and theresponse matrix is found by dividing all of these distribu-tions by the number of incident photons in the respectivecause bin.The absorption of photons in the atmosphere dependson the zenith angle of the observation and the altitude ofthe balloon, both of which vary with time. Therefore, anobservation of a source has to be considered as a sequenceof observations I = 1 , . . . , n obs of duration T I , describedby individual response matrices R I . Assuming that thesource spectrum N cj is constant during the entire obser-vation time T = (cid:80) I T I , the detector response can be de-scribed by an average response matrix: N ei = n c n cφ (cid:88) j =1 (cid:32) T n obs (cid:88) I =1 T I R ij,I (cid:33) N cj . (10)Since the individual observations only differ by the absorp-tion in the atmosphere, the response matrices will onlydiffer by a factor in each column: R I = R diag( t I ) . (11)Therefore, absorption in the atmosphere as a function ofphoton energy in time interval I , t I , can be separated fromthe other effects listed above. The vector t I contains theprobability that a photon in the j -th true energy bin willnot be absorbed in the atmosphere: t I,j = exp( − d µ ( E j ) /ρ ) , (12)where d is the slant depth in g cm − and µ ( E j ) /ρ is themass attenuation coefficient in air within energy bin E j .Thus, it is sufficient to simulate the detector response ex-cluding atmospheric absorption once ( R ) and then applythe absorption effects for each time interval. The optimal number of unfolding iterations has to bedetermined before attempting to unfold the experimentaldata in order to avoid potential biases. Generally, the un-folded spectrum becomes a better representation of thetrue spectrum with increasing number of iterations. How-ever, at the same time fluctuations of the data are am-plified leading to the same problem introduced by simplyinverting the response matrix. This can be understoodsince in some way the unfolding matrix from Eq. (7) is an“inverse” of the response matrix. Since the response matrixsmears fluctuations due to the limited detector response,the unfolding matrix will have the opposite effect, ampli-fying fluctuations. In addition each new prior is basedon the unfolding result from the previous iteration, whosefluctuations in this way feed back into the unfolding loop.The iteration should thus be terminated when a compro-mise between a good representation of the true spectrumand small artificial fluctuations is reached.One way to achieve this is to fold a hypothetical in-put spectrum and azimuthal distribution with the responsematrix and then sample the same number of events as inthe experimental data set from this distribution. These“ measured ” data are then unfolded to find the number ofiterations that yields the best agreement between the re-sult and the known true input. This should be repeated forvarious different inputs to find an optimum independentof the true spectrum and polarization fraction or angle.The disadvantage of this method is that the ideal num-ber of iterations may depend considerably on whether thetrue spectrum is smooth or contains features such as emis-sion or absorption lines. The choice of number of iterationsmay therefore lead to a significant systematic bias.A better approach is to use the convergence of the un-folding process itself to determine when to stop the it-eration. After each iteration, the unfolded spectrum isfolded with the response matrix, ˜ N e ( r ) i = (cid:80) j R ij ˆ N c ( r ) j ,and compared to the measured data, N ei . A convergencecriterion can then be defined using the change in χ be-tween ˜ N e ( r ) i and the measured data between two itera-tions r and r + 1 [15]: ∆ ˜ χ ( r, r + 1) = χ ( ˜ N e ( r ) , N e ) − χ ( ˜ N e ( r +1) , N e ) . (13)This quantity decreases monotonically during the itera-tion process. However, after a large number of iterationsas ∆ χ ( r, r + 1) → the unfolding would introduce thesame large-scale fluctuations as simply inverting R i,j . Toavoid this, the iteration is terminated once ∆ χ ( r, r + 1) falls below a certain value ∆ χ . The value of this limitcan be determined beforehand using a simple Monte Carlosimulation similar to what would be used to determine afixed maximum number of iterations as described above.The best value of ∆ χ generally depends less on theshape of the true spectrum than the number of iterations.4 .4. Statistical uncertainties of the unfolded spectrum Statistical uncertainties of the measured data are prop-agated through the unfolding process using a bootstrapmethod: the number of events within each bin of the in-put data is varied randomly r = 1 , . . . , n times accordingto a Poisson distribution with the measured value as themean, and the unfolding is repeated. The statistical errorin bin j is then determined by comparing each unfoldingresult ˆ N c ( r ) j to the average (cid:104) ˆ N c (cid:105) j : ( σ cj ) = 1 n − n (cid:88) r =1 (cid:0) ˆ N c ( r ) j − (cid:104) ˆ N c (cid:105) j (cid:1) . (14)Likewise, bin-to-bin correlations are obtained: cov( i, j ) = 1 n n (cid:88) r =1 (cid:0) ˆ N c ( r ) i − (cid:104) ˆ N c (cid:105) i (cid:1)(cid:0) ˆ N c ( r ) j − (cid:104) ˆ N c (cid:105) j (cid:1) . (15) For the studies in the next sections we made the as-sumption that the azimuth angle of the scattered photonscan be determined perfectly. This means that the responsematrix is diagonal in azimuth, which makes it a sparse ma-trix, essentially reducing its size by one dimension. Theauthors have, therefore, implemented the unfolding algo-rithm using sparse matrix routines from the
SuiteSparse package [16], which lead to a significant performance im-provement, compared to an implementation using densematrix algorithms. Even when taking into account uncer-tainties in the azimuth distribution, the response matrixwill very likely still be sparse, meaning sparse matrix al-gorithms are likely a good choice in any case.
4. Application to simulated X-Calibur data
In this section we will apply the method described be-fore to a simulation of a -hour observation of the Crabnebula with the polarimeter X-Calibur flown in the focalplane of the InFOC µ S telescope.
The detection principle of X-Calibur is illustrated inFig. 2. Incoming X-rays Compton-scatter in a
14 cm long,
13 mm diameter scintillator rod, and the scattered photonsare detected by 32 and thick Cadmium-Zinc-Telluride (CZT) detectors, each occupying an area of ×
20 mm divided into × pixels [5]. In order to reducesystematic uncertainties, X-Calibur continuously rotatesaround the optical axis. Observables are the azimuthalangle φ and energy E of the scattered photon, and thecoordinate z along the optical axis at which the scatteredphoton was detected.The z -coordinate of the detected photon is related tothe scattering angle θ and the depth along the axis ofthe scatterer at which the Compton scattering occurred,which, however, cannot be measured in the experiment. PMT
16 cmCZT (5 & 2 mm thick)Scintillator ( ⌀
13 mm) X-rays0 z Figure 2: Principle of the X-Calibur polarimeter. A grazing inci-dence X-ray mirror focuses the X-rays onto the polarimeter. In thissketch the X-rays enter the detector from the right and scatter off ascintillator slab. The scattered X-rays are then detected by one ofthe CZT detectors that surround the scintillator. A coincidence re-quirement between the trigger from the photomultiplier tube at theend of the scintillator and events in the CZTs can be used to suppressbackgrounds, but is not used in the analysis presented here.
Therefore, z carries important information on the loss ofenergy in the scatterer. Including this variable in the anal-ysis, therefore, provides valuable information helping toreconstruct the energy of the incident photon.In order to test the unfolding method, the response ofthe X-Calibur polarimeter to incoming X-rays was sim-ulated using Geant4 [17, 18]. The response of the CZTdetector was then simulated by propagating a charge pro-portional to the amount of energy deposited inside thedetector through the electric field of the detector and inte-grating the induced charge at the electrodes. Additionally,the electronic response was smeared with a Gaussian dis-tribution with a resolution of . The detector thresholdwas modeled as a hard cutoff at
20 keV measured energy.The X-ray mirror was simulated by assigning a weight toeach event proportional to the effective collection area atthe incident photon energy.Three data sets were generated. The first (
DS1 ) con-sisted of horizontally polarized photons with an E − spectrum in the range between and
85 keV to be usedto extract the response matrix. It took into account thedetector response and mirror effective area as describedabove but no atmospheric absorption.The second dataset (
DS2 ) was used to test the unfold-ing. It consisted of the same detector simulation but witha steeper incident spectrum, E − . , as found in the CrabNebula [19]. Additionally, observation time dependent ab-sorption of photons in the atmosphere (due to the changeof source elevation) was introduced to simulate the obser-vation of an astrophysical source (in our example the CrabNebula as observed from Ft. Sumner, NM, USA) at an at-mospheric overburden of . − . The absorption as afunction of energy was calculated from a spline interpo-lation of the NIST XCOM tables for air [12]. The polar-ization of the Crab Nebula has been measured at . , . [2] and >
100 keV [20, 21]. The polarization seemsto increase from ∼ at . and . to values of [20] or even > [21] at >
100 keV . We assume an5 ime of day / hrs2 3 4 5 6 7 8 9 10 11 R a t e / H z ° E l e v a t i on / Figure 3: Simulated event rate in X-Calibur (thin line) and elevationof the Crab pulsar (bold line) as a function of time on 9/24/2014. energy dependent polarization fraction p ( E ) of: p ( E ) = 0 . − . (cid:16) E/ keV − (cid:17) + 1 , (16)which leads to p ∼ at
20 keV and p ∼ above
70 keV with a smooth transition around
40 keV . The po-larization angle does not change in this example. Thenumber of events was chosen to match an -hour observa-tion of the Crab Nebula with X-Calibur assuming a fluxof F ( E ) = 10 . (cid:18) E (cid:19) − . cm − s − keV − (17)as measured by the Swift Burst Alert Telescope [19]. Eachevent in this dataset was assigned a time of day, whichprogressed to match the expected rate. The zenith an-gle of the source at this time was stored with each event.The observation time dependence of the event rate that re-sulted from this zenith angle variation is shown in Fig. 3.This dataset contained a total of 23,872 events.A third dataset ( DS3 ) was created from the second oneto include an absorption line at an energy of E line = 40 keV by discarding events with a probability of p absorb ( E phot ) = 0 . × exp (cid:18) − ( E phot − E line ) σ (cid:19) (18)with a width of σ line = 5 keV , which is comparable tothe resolution of the simulated instrument. Figure 4shows the spectra of measured photon energies E det fordatasets DS2 and
DS3 as well as their ratio, which clearlyshows the effect of the absorption line.
From data set
DS1 the response matrix R was con-structed as described in Section 3.2. Using the tables in E v en t s
10 Power lawPower law + absorption line / keV det
E20 30 40 50 60 70 80 90 R a t i o Figure 4: Spectrum of measured photon energies in data sets
DS2 (pure power-law) and
DS3 (power-law + absorption line). The bot-tom panel shows the ratio of the two spectra highlighting the effectof the absorption line. The large bin-to-bin variations at low energiesare a binning effect: the discretization of energies in the detector islinear whereas the binning in these histograms was done on a loga-rithmic scale.
Ref. [12] and the elevation angle of the Crab (see Fig. 3),the atmospheric absorption in each energy bin was calcu-lated in intervals. These absorption coefficients werethen applied to the response matrix R using Eqns. (10)and (11).In this example we consider the simplified case wherethe response is diagonal in azimuth, i. e. we assume φ e = φ c , and measured and “true” azimuth angles will be binnedin the same way, n cφ = n eφ =: n φ . The azimuth needs to beincluded in the response matrix in order to reconstruct thedependence of the azimuth distribution on the true energyin the unfolding. Assuming φ e = φ c is not a conceptualrestriction and can be removed easily. It is merely a simpli-fication used in the study of X-Calibur data in this paperbecause the precision at which the scattering angle is mea-sured does not limit the reconstruction of the polarizationangle.A representation of the energy response of X-Caliburis shown in Fig. 5. On average, the detected energy is lessthan the incident photon energy as the photon transferssome of its energy to the Compton electron. The distri-bution of measured energies shows a long tail towards lowvalues. This is due to the combined effect of the physicsof Compton scattering, fluorescence photons escaping theCZTs, the depth-dependent charge collection efficiency ofthe CZT detectors, and occasional Compton scatteringrather than photoelectric effect interactions in the CZTcrystals. At high energies, the efficiency is dominated bythe small effective area of the mirror, while at low energies6 / keV ph E20 30 40 50 60 70 / k e V de t E P r obab ili t y Figure 5: Projection of the response matrix onto the E detected -vs.- E photon plane. absorption in the atmosphere is the dominant effect.Finally, the modulation factor was determined as afunction of energy. In each energy bin the azimuthal scat-tering angle distribution measured for a polarizedbeam was fitted with a sine function, f ( χ ) = A sin( nχ − φ ) + B, (19)with the free parameters: amplitude A , periodicity π/n ,phase φ , and offset B . The modulation factor is then µ = A/B. (20)No energy dependence of µ was found and an averagevalue of (cid:104) µ (cid:105) = 0 . ± . was determined. In order to determine the optimal number of iterations,the unfolding was applied to simulated data distributions.A true energy spectrum and energy-dependent scatteringangle distribution was folded with the response matrix de-termined above. The result is the 3-dimensional distri-bution of detected energies, scattering angles and z coor-dinates one would find in the absence of statistical errors.This distribution was then scaled such that its integral cor-responded to the total number of events in the test datasetto be unfolded ( DS2 and
DS3 ).The event number N in each bin was then random-ized according to a normal distribution with mean N andstandard deviation √ N . This simulated dataset was thenunfolded with the same prior that was going to be used tounfold the simulated data. After each iteration• the statistical errors of the unfolded spectrum werecalculated according to Eq. (14);• ∆ ˜ χ was calculated according to Eq. (13);• and the χ between the unfolded and the true spec-trum was calculated. (r+1, r) χ∼∆ χ
10 spectrum True E spectrum True E spectrum True EStop condition
Figure 6: Difference between true and unfolded spectrum as a func-tion of ∆˜ χ , defined in Eq. (13), as described in Section 3.3. Theunfolding proceeds from right to left as ∆˜ χ decreases monotoni-cally. At first, χ decreases until reaching a minimum after whichit increases again due to numerical noise. The different colors repre-sent true spectra with three different powerlaw indices. In all casesan energy-independent polarization fraction of . was assumed. Theprior was in all cases an E − spectrum flat in azimuth. The dottedlines represent individual unfolded spectra. For each true spectrum30 randomizations are shown. The thick lines are the average overall 100 variations per true spectrum. When applied to data the iter-ation will proceed until ∆˜ χ drops below the value indicated by thevertical dash-dotted line labeled “Stop condition”. This process was repeated times.The result of this procedure for three different truespectra with an energy-independet polarization fractionof . is shown in Fig. 6. In all cases an E − prior witha flat azimuthal distribution was chosen. Generally, theagreement between the unfolded and the true spectrumimproves with increasing number of iterations, i. e. decreas-ing ∆ ˜ χ . However, when exceeding a certain number ofiterations, the χ starts to vary by many orders of mag-nitude. At this point, the unfolding procedure starts toamplify the statistical fluctuations. Generally, the unfold-ing iteration should be stopped before reaching this region.More detailed analysis shows that the greatest improve-ment is reached within the first few iterations.Two unfolding tests were performed, both based on thesame input described in Sect. 4.1. First, the data wereunfolded starting from an E − prior spectrum, then withan E − prior. Assuming the true spectrum was unknown,the optimal number of iterations for three different truespectra each with polarization fractions of , . , and was determined using the method described in this section.Furthermore, spectra with an absorption line with a widthof at E = 40 keV were studied. The decision on thenumber of iterations for the two examples was then reachedbased on curves similar to those in Fig. 6. We chose toiterate until ∆ ˜ χ < , which resulted in the number ofiterations listed in Table 1 for the four cases studied here.7 able 1: Number of unfolding iterations resulting from the stoppingcondition ∆˜ χ < for the two priors and the two true spectra. Prior index: − − Power-law
Power-law + line
To test the unfolding method, the simulated data sets
DS2 and
DS3 were first unfolded starting from a prior thatfollowed an E − energy spectrum and was flat in azimuthfor all energies. In a second test, an E − prior, againflat in azimuth, was used. Figure 7 shows the spectra ofmeasured energies, the true spectra and the energy spectraobtained as a result of the unfolding.Generally, the energy spectrum is reproduced reason-ably well. The lowest and highest energy bins are problem-atic due to the low detection efficiencies (see the plots inthe right column of Fig. 7). The absorption line at
40 keV is smeared out visibly, as seen in Figures 7(b) and (d).However, one should note that the width of the line isvery close to the energy resolution of the simulated CZTdetectors and there is additional energy smearing due toenergy losses in the Compton scattering, which cannot bemeasured in the experiment. It is important to note thatin both cases (pure power-law and power-law plus absorp-tion line) the choice of prior does not have a significantinfluence on the result beyond the statistical fluctuations.In a separate step after the unfolding the azimuthaldistribution in each energy bin was fitted with a sine func-tion, Eq. (19), in order to determine the polarization frac-tion and angle as a function of energy. In order to improvethe accuracy of the determination of polarization param-eters, energy bins were combined into larger bins. Theerror bars were determined from the uncertainties of thefit parameters through standard error propagation. As anexample, the azimuthal distributions and fits for the caseof a power-law with absorption line, starting from an E − prior are shown in Fig. 8.The resulting polarization fraction as a function of en-ergy is shown in Fig. 9 for all four cases (pure power-law spectrum, power-law plus absorption line, each withan E − and a E − prior). In most cases, the change inpolarization fraction would be observed by a measurementsuch as that simulated for this study. At the highest ener-gies the statistical uncertainties become rather large dueto the low statistics in each energy bin. The results in theindividual energy bins are somewhat positively correlateddue to the mixing of events in neighboring bins throughthe unfolding, which leads to slightly larger error bars thanone would expect from the scattering of points. Further-more, it should be noted that a χ test reveals a slightlybetter agreement of the results with the model in case ofthe E − prior ( χ /N df = 11 . / versus . / and . / versus . / ). Azimuth / degrees0 50 100 150 200 250 300 350 deg s c m P ho t on s / k e V
10 19.3keV22.8keV26.9keV31.7keV35.8keV38.9keV42.2keV47.8keV58.7keV
Figure 8: Reconstructed azimuthal distributions in eight energy binswith sine fits for the case of a power-law with absorption line spec-trum. The unfolding was started from an E − prior. The numberson the right indicate the central energy of the energy bins. E / keV30 40 50 60 70 80 P o l a r i z a t i on F r a c t i on Power law, E prior Power law, E prior Power law + absorption line, E prior Power law + absorption line, E
Figure 9: Reconstructed polarization fraction as a function of pri-mary energy after the unfolding prodecure starting from an E − andan E − prior spectrum model polarization fraction, Eq. (16). s c m P ho t on s / k e V
10 TrueRecontructedE / keV20 30 40 50 60 70 R a t i o E v en t s
10 DataRefolded / keV det
E20 30 40 50 60 70 80 90 R a t i o (a) Power-law starting from E − prior. s c m P ho t on s / k e V
10 TrueRecontructedE / keV20 30 40 50 60 70 R a t i o E v en t s
10 DataRefolded / keV det
E20 30 40 50 60 70 80 90 R a t i o (b) Power-law with absorption line starting from E − prior.Figure 7: Result of the unfolding prodecure starting from an E − prior spectrum for (a) the case of a pure power-law spectrum and (b) apower-law with an absorption line at an energy of
40 keV . Figures (c) and (d) show the same but starting from an E − prior spectrum. Left:
Unfolded energy spectrum compared with the true spectrum (see Eq. (17) for the power-law component).
Right:
As a cross-check of theunfolding procedure, the results were folded with the response matrix and compared to the input data. The figures show a projection of thedata onto the axis of measured energy. s c m P ho t on s / k e V
10 TrueRecontructedE / keV20 30 40 50 60 70 R a t i o E v en t s
10 DataRefolded / keV det
E20 30 40 50 60 70 80 90 R a t i o (c) Power-law starting from E − prior. s c m P ho t on s / k e V
10 TrueRecontructedE / keV20 30 40 50 60 70 R a t i o E v en t s
10 DataRefolded / keV det
E20 30 40 50 60 70 80 90 R a t i o (d) Power-law with absorption line starting from E − prior.Figure 7: See caption on page 9 for details. . Summary Experiments like X-Calibur, PoGOLite, GEMS, andXIPE are able to measure polarization of X-rays from as-trophysical sources over a wide range of energies and en-able spectropolarimetric measurements. Typically, resolu-tion, effective area, and energy response of X-ray polarime-ters depend strongly on the energy of the incident X-ray.Due to the generally steep source spectra, these energydependencies need to be taken into account when recon-structing the flux, polarization fraction and polarizationdirection energy spectra.There are several methods that can be used to achievethis. Generally, the distributions of input and output vari-ables are related through a response matrix, which de-scribes the probability distribution of outputs for each setof input parameters. Forward folding is a technique com-monly used if a model of the true spectrum should befitted to the data. The model spectrum is folded withthe response matrix and compared to the measured data,and model parameters are adjusted to achieve a best fit.The disadvantage of this method is that it cannot delivermodel-independent data points. By combining forwardfolding with Markov Chain Monte Carlo, the contents ofeach bin can be treated as an independent parameter al-lowing a model-independent spectral reconstruction. How-ever, this becomes very computationally expensive if thenumber of bins is large, e. g. in case of multi-dimensionalproblems. The inverse of the response matrix can be usedto transform the measured distributions into estimates ofthe true spectrum. However, direct matrix inversion usu-ally leads to an amplification of statistical fluctuations andconsequently a wildly fluctuating result where neighboringvalues have a strong negative correlation.In this paper an unfolding-based method is describedthat makes it possible to reconstruct the flux energy spec-trum and to infer the polarization fraction and polariza-tion direction as function of the true energy of the de-tected photons. The method makes use of a Bayesianmulti-dimensional iterative unfolding procedure [14]. Fora scattering or photo-effect polarimeter, the input param-eters typically are the distributions of measured photonenergy and azimuthal scattering angle or azimuth angleof the emitted photo-electron. The output usually is atwo-dimensional distribution of true photon energies andscattering angles. However, additional parameters can eas-ily be included both on the input and the unfolded data.In Section 4, this is demonstrated for an additional inputparameter.Starting from a prior distribution, which can reflect thebest knowledge about or a best guess of the true spectrumand polarization, the data are then unfolded iteratively asdescribed in Section 3. Statistical uncertainties and thecovariance matrix of the unfolding result are computed byvarying the input within its statistical errors, and observ-ing variations and correlations of the output spectrum. InSection 4 this method was applied to a simulated dataset based on the X-Calibur scattering polarimeter. Two sce-narios were studied: first, the case of a simple power-lawspectrum; secondly, a Gaussian shaped absorption line wasadded at an energy of
40 keV . In both cases the true spec-trum was reconstructed reasonably well within 6 or lessiterations, in particular when considering that the widthof the absorption line is close to the resolution of the simu-lated detector. After the unfolding the azimuthal scatter-ing angle distributions were fitted with sine functions inbins of “true” energy yielding the polarization fraction andangle as a function of energy. It was shown that the resultdoes not depend strongly on the chosen prior spectrum,meaning the choice of prior does not lead to a strong biasof the result.
Acknowledgements
The authors are grateful for NASA funding from grantsNNX10AJ56G & NNX12AD51G as well as discretionaryfunding from the McDonnell Center for the Space Sciences.
ReferencesReferences [1] H. Krawczynski et al., Scientific Prospects for Hard X-ray Po-larimetry, Astropart. Phys. 34 (2011) 550–567.[2] M. C. Weisskopf et al., A Precision Measurement of the X-RayPolarization of the Crab Nebula without Pulsar Contamination,Astrophys. J. 220 (1978) L117–L121.[3] J. E. Hill et al., The design and qualification of the GEMS x-ray polarimeters, in: Proc. SPIE Space Telesc. Instrum. 2012:Ultrav. to Gamma Ray, Vol. 8442, 2012.[4] P. Soffitta et al., XIPE: the X-ray imaging polarimetry explorer,Exp. Astron. 36 (2013) 523–567.[5] Q. Guo et al., Optimization of the design of the hard X-raypolarimeter X-Calibur, Astropart. Phys. 41 (2013) 63–72.[6] M. Pearce et al., Balloon-borne hard X-ray polarimetry withPoGOLite, Proc. 2012 IEEE NSS. arXiv: astro-ph/1211.5094.[7] W. Heitler, The Quantum Theory of Radiation, ClarendonPress, Oxford, 1936.[8] R. D. Evans, The atomic nucleus, McGraw-Hill, New York,1955.[9] J. Sánchez Almeida, V. Martínez Pillet, Polarizing propertiesof grazing-incidence x-ray mirrors: comment, Appl. Optics 32(1993) 4231–4235.[10] J. Katsuta et al., Evaluation of polarization characteristics ofmultilayer mirror for hard X-ray observation of astrophysicalobjects, Nucl. Instrum. Meth. Phys. Res. A 603 (2009) 393–400.[11] J. Tueller et al., InFOCµS Hard X-ray Imaging Telescope, Exp.Astron. 20 (2005) 121–129.[12] J. H. Hubbell, S. M. Seltzer, Tables of X-Ray Mass AttenuationCoefficients and Mass Energy-Absorption Coefficients, Online(version 1.4) 2004, Retrieved 10/18/2013. Originally publishedas NISTIR 5632 (1995).URL http://physics.nist.gov/xaamdi [13] P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems,SIAM, Philadelphia, 1998.[14] G. D’Agostini, A multidimensional unfolding method based onBayes’ theorem, Nucl. Instrum. Meth. Phys. Res. A 362 (1995)487–498.[15] H. Ulrich, Untersuchungen zum primären Energiespektrum derkosmischen Strahlung im PeV-Bereich mit dem KASCADE-Experiment, Ph.D. thesis, Universität Karlsruhe, in German(2004). RL [16] Y. Chen, T. A. Davis, W. W. Hager, S. Rajamanickam,Algorithm 8xx: CHOLMOD, supernodal sparse Choleskyfactorization and update/downdate, Tech. Rep. TR-2006-05,CISE Dept., Univ. of Florida, Gainesville, FL (2006).URL [17] S. Agostinelli et al., Geant4 — a simulation toolkit, Nucl. In-strum. Meth. A 506 (2003) 250–303.[18] J. Allison et al., Geant4 developments and applications, IEEETrans. Nucl. Sci. 53 (2006) 270–278.[19] J. Tueller et al., The 22 Month Swift-BAT All-Sky Hard X-RaySurvey, Astrophys. J. Suppl. S. 186 (2010) 378–405.[20] A. J. Dean et al., Polarized Gamma-Ray Emission from theCrab, Science 321 (2008) 1183–1185.[21] M. Forot et al., Polarization of the Crab Pulsar and Nebulaas Observed by the INTEGRAL/IBIS Telescope, Astrophys. J.688 (2008) L29–L32.[17] S. Agostinelli et al., Geant4 — a simulation toolkit, Nucl. In-strum. Meth. A 506 (2003) 250–303.[18] J. Allison et al., Geant4 developments and applications, IEEETrans. Nucl. Sci. 53 (2006) 270–278.[19] J. Tueller et al., The 22 Month Swift-BAT All-Sky Hard X-RaySurvey, Astrophys. J. Suppl. S. 186 (2010) 378–405.[20] A. J. Dean et al., Polarized Gamma-Ray Emission from theCrab, Science 321 (2008) 1183–1185.[21] M. Forot et al., Polarization of the Crab Pulsar and Nebulaas Observed by the INTEGRAL/IBIS Telescope, Astrophys. J.688 (2008) L29–L32.