Pulsar Emission above the Spectral Break - A Stacked Approach

Pulsar Emission above the Spectral Break - A Stacked Approach


NASA’s Fermi space telescope has provided us with a bountiful new population of gamma-ray sources following its discovery of over 150 new gamma-ray pulsars. One common feature exhibited by all of these pulsars is the form of their spectral energy distribution, which can be described by a power law followed by a spectral break occurring between 1 and 8 GeV. The common wisdom is that the break is followed by an exponential cutoff driven by radiation-reaction-limited curvature emission. The discovery of pulsed gamma rays from the Crab pulsar, the only pulsar so far detected at very high energies (E100 GeV), contradicts this “cutoff” picture. Here we present a new stacked analysis with an average of 4.2 years of data on 115 pulsars published in the 2nd Fermi-LAT catalog of pulsars. This analysis is sensitive to low-level 100 GeV emission which cannot be resolved in individual pulsars but can be detected from an ensemble.

I Introduction


Figure 1: SEDs data points and flux upper limits for the Geminga pulsar. Measurements of the Crab Nebula and pulsar are plotted for comparison. It is clear, even in the phase-resolved analysis, that the SED falls slower than an exponential and appears more consistent with a simple power-law. Figure taken from Aliu et al. (2015).

One common feature exhibited by all known gamma-ray pulsars is the form of their spectral energy distribution (SED) which can be described by a power-law followed by a spectral break occurring between 1 and 8 GeV (Abdo et al., 2013). The unanimity of the break energy across the entire Fermi-LAT pulsar sample is suggestive that the sites of acceleration and processes of gamma-ray emission are common across different pulsar types and that they are not strongly dependent on the pulsar spin or energetics. Further, it has been shown that across the Fermi-LAT pulsar sample the spectral-break energy is weakly correlated with the magnetic-field strength at the light cylinder (Abdo et al., 2010b, 2013). Such behavior is expected in models where emission is produced by curvature radiation (CR) occurring at the radiation-reaction limit in the outer magnetosphere (Harding et al., 2008; Abdo et al., 2010b). This has become the most favored general description of gamma-ray emission from pulsars in the Fermi-LAT era. In these models one expects that the SED will fall off exponentially above the break energy. There is, however, compelling evidence suggesting that CR occurring in the outer magnetosphere is not a complete description of pulsar emission at, and above, the GeV SED break:

  1. The discovery of power-law-type emission from the Crab pulsar at energies exceeding 100 GeV1 cannot be easily explained by curvature emission from the outer magnetosphere (Aliu et al., 2011; Lyutikov et al., 2012) unless the radius of curvature of the magnetic field line is larger than the radius of the light cylinder (Bednarek, 2012). Some recent models attribute the pulsed very-high-energy (VHE; 100 GeV) emission from the Crab pulsar to inverse-Compton (IC) scattering originating in the outer magnetosphere (Lyutikov et al., 2012; Du et al., 2012; Lyutikov, 2012) or to IC scattering from beyond the light cylinder (Aharonian et al., 2012; Pétri, 2012).

  2. The radiation-reaction limit of CR occurs when the acceleration gains achieved by an electron are equaled by radiation losses. The photon energy at which this occurs in the outer magnetosphere can be expressed in terms of the pulsar period, the surface magnetic field strength, the radius of curvature of the accelerated particle and an efficiency factor. Lyutikov et al. (2012) has shown that the break-energy values for several pulsars reported in the first Fermi-LAT pulsar catalog are so high that they require the efficiency factor and radius of curvature to approach or even reach their maximal allowable values2.

  3. Recent studies of the Geminga pulsar with Fermi-LAT and VERITAS (see Figure 1) show that the SED above the GeV break is compatible with a steep power law (Lyutikov, 2012; Aliu et al., 2015), but no emission has been seen above 100 GeV. Similar conclusions can be drawn from an analysis of the Vela pulsar with Fermi-LAT data from Leung et al. (2014), who show that multi-zone or time-dependent emission models are needed to fit the slower-than-exponential fall of the SED above 10 GeV.

The question of whether the Crab pulsar is unique, or whether non-exponentially-suppressed gamma-ray spectra are common in gamma-ray pulsars is of great importance. Beyond the modeling of pulsar emission, questions concerning the emission spectra of pulsars have significant implications for galactic dark matter searches, where unassociated gamma-ray excesses can be interpreted as the remnants of dark matter annihilation (e.g., Abazajian & Kaplinghat 2012). Since pulsars are likely the main background for these searches, categorizing the shape of pulsar spectra is a critical step towards validating any indirect dark matter signal in the gamma-ray domain. To search for non-exponentially-suppressed emission above 50 GeV, we have performed a stacked analysis of gamma-ray pulsars which is sensitive to emission which cannot be resolved in the Fermi-LAT analysis of individual objects, but can be detected if aggregated from an ensemble. A stacked analysis which yields evidence of cumulative emission above 50 GeV would prove that some population of gamma-ray pulsars clearly exhibits non-exponentially-suppressed emission. This would indicate that inverse-Compton or wind-zone emission is common in gamma-ray pulsars and that pulsars contribute to the sub-TeV diffuse emission of the galaxy.

Ii Analysis

Figure 2: Aperture photometry analysis steps for the Crab pulsar. Panel (a) plots the phase distribution (light curve) of the Crab pulsar from 5.2 years of Fermi-LAT observations. The Off phase range, [0.71 0.99], is defined in the 2nd Fermi-LAT catalog of gamma-ray pulsars (2PC). Panel (b) plots the distribution of photon energies for events which fell in the On and Off phase ranges. The Off events have been scaled by which is the ratio of the On phase gate(s) size to the Off gate(s) size . Panel (c) shows the energy distribution of the excess events and panel (d) shows the significance of the excess in each energy bin. Panel (e) shows the Fermi-LAT exposure for the ROI used in each energy bin determined from gtexposure. In panel (f) the Crab pulsar AP SED is plotted alongside the Crab pulsar SED determined from a likelihood fit done in the 2PC. A broken power-law fit to Fermi-LAT and VERITAS data from Aliu et al. (2011) is plotted, as well as the VERITAS 100 GeV bow-tie. Below the SED plotted in panel (f) is the ratio of the AP flux to the 2PC flux in each bin, showing the level of agreement between the AP method and the likelihood method. This figure is taken from McCann (2014).

ii.1 The Aperture Photometry Method

A maximum likelihood fitting procedure is typically employed when performing spectral analysis of Fermi-LAT data. The Fermi-LAT data can also be analyzed with an Aperture Photometry (AP) method where the raw event counts from a region of interest (ROI) are combined with a measure of the instrument exposure (cm s) to the region to determine the flux. This AP method is less sensitive and less accurate than the likelihood fitting procedure but it “provides a model independent measure of the flux” and it “is less computationally demanding”3. We demonstrate here that the AP method can be used to produce accurate SEDs from multi-year pulsar data sets since an accurate determination of the background rate can be measured from the “Off-pulse” phase range. The analysis presented here proceeds as follows:

  1. Over the 100 MeV to 1 TeV energy range, logarithmically-spaced energy binning with 4 bins per decade is chosen.

  2. An ROI is chosen around each pulsar with an energy-dependent radius. The radius chosen is three times the 68% point-spread-function (PSF) containment radius determined from a “front-conversion” Vela analysis by Ackermann et al. (2013). In order to maintain sufficient statistics at high energies, the radius of the ROI was fixed to 0.45 above 10 GeV.

  3. The Fermi-LAT analysis tools gtselect, gtmktime, gtbin and gtexposure are then run over each pulsar ROI for all observations performed within the period of validity of the pulsar timing solution.

  4. The photon event list is barycentered and phase-folded using the Tempo2 package (Hobbs et al., 2006) with the Fermi Tempo2 plugin and the corresponding timing solution.

  5. Within each energy bin, a cut on phase is applied and events which fall within the Off phase region and those which fall outside this region - the On phase region - are selected. The ratio of the size of the On phase range to the size of the Off phase range, defined as , is then used to scale the number of event counts in the Off phase region (N) to the number in the On region (N).

  6. The number of excess pulsed events is then defined as NNN and the flux is N divided by the exposure () calculated in step 3 using gtexposure. The significance of the excess is calculated using Equation 17 from Li & Ma (1983).

Following this procedure one can derive the energy distributions for the On and Off phase regions, and the instrument exposure, for any pulsar. These distributions and the derived AP SED are shown for the Crab pulsar in Figure 2.

ii.2 Stacking The Pulsar Data Sets

Fermi-LAT has detected over 150 new gamma-ray pulsars4 and the stacking performed in this work uses 115 pulsars listed in the Second Fermi-LAT Catalog of Gamma-ray Pulsars (Abdo et al., 2013), which shall be referred to as 2PC throughout5. The 115 pulsar sample is composed of 39 millisecond pulsars and 76 “young” non-recycled pulsars with an average data set spanning 4.2 yr6. The six AP analysis steps listed in Section II.1 were followed for each pulsar and using the resulting values of N, N and , it is quite simple to determine the total excess,


the total exposure,


and thus, the average flux,


for N pulsars in a given energy bin, . The significance of the total excess is determined by the generalized version of Equation 17 from Li & Ma (1983) (see Aharonian et al. 2004). In cases where the significance is less than 2, the method of Helene (1983) is used to derive the 95% confidence-level upper limit on the total excess, which is in turn used to compute a flux upper limit.

Iii Results

Figure 3: Panel (a) shows the average flux (square markers) from 76 young pulsars determined by dividing the total excess by the total exposure (see Equations 13). The dashed-line histogram shows one over the total exposure, indicating the flux which would correspond to a single excess photon. This is the minimum possible flux which could be measured given the total exposure. The gray cross shows the most constraining limit on emission from a single pulsar in the 56.2100 GeV range presented in the 2PC. The 2PC presented no limits at higher energies. The broken power-law fit to the Crab pulsar data from Aliu et al. (2011) is plotted for scale. Panel (b) plots the same quantities for the stacked analysis of 39 millisecond pulsars. This figure is adapted from figures presented in McCann (2014).

The stacking analysis results for the young pulsar and millisecond pulsar ensembles are shown in Figure 3. No significant excesses are seen in these analyses at energies above 50 GeV. Upper limits on the average flux, determined at the 95% confidence-level, are listed in Table 1 for three energy bins above 50 GeV. Limits are also presented in units of the Crab pulsar where the broken power-law fit to the Fermi-LAT and VERITAS data presented in Aliu et al. (2011) defines a Crab pulsar unit. In addition to these analyses, we stacked sub-samples of the data where each sub-sample was composed of the 10 pulsars with the largest value of a given parameter. Sub-sample selections based on gamma-ray luminosity, spin-down power, spin-down power over distance squared, gamma-ray photon flux and non-thermal X-ray energy flux were investigated7. No significant excesses were observed above 50 GeV in any of these sub-sample stacking analyses.

The shape of the average young pulsar and average millisecond pulsar SEDs were categorized by fitting a power law times a super-exponential cutoff function {linenomath}


to the SED data. These fits are presented in Figure 4. Fixing reduces Equation 4 to a power law times an exponential cutoff function and, as expected, this functional form does not reproduce the sub-exponential fall of the SED above the break. However it can be used to measure the average flux-weighted value of the spectral index () and cutoff () parameters (Abdo et al., 2013). It is clear from Figure 4 that the average SEDs have qualitatively the same shape, with the average flux from the 39 millisecond pulsars about an order of magnitude lower than the average flux from the 76 young pulsars. The spectral parameters derived from fitting the two ensembles are remarkably similar and are presented in the caption of Figure 4.

All Young Pulsars Millisecond Pulsars
Energy Range Flux Limit Flux Limit Flux Limit Flux Limit Flux Limit Flux Limit
[GeV] [ [Crab pulsar [ [Crab pulsar [ [Crab pulsar
cms] units] cms] units] cms] units]
56.2 100 1.57 0.07 2.03 0.09 1.44 0.07
100 177 1.52 0.31 1.88 0.38 1.14 0.23
177 316 1.34 1.21 1.96 1.76 0.50 0.45
Table 1: Limits at the 95% confidence level on the average flux from stacked ensembles of gamma-ray pulsars. The limit values presented in Crab pulsar units assume the broken power-law fit to the Crab pulsar data from Aliu et al. (2011) is a Crab pulsar flux unit.
Figure 4: The average SEDs derived from the stacking of the 76 young pulsars and 39 millisecond pulsars. The SEDs are each fit with a power law times a super-exponential cutoff keeping both fixed to unity and allowing it to float. For the pure exponential cutoff case () the best fit value is 0.540.05 for the millisecond pulsars and 0.410.01 for the young pulsars while the best fit values are 3.600.21 GeV and 3.540.04 GeV, respectively. Allowing to float we find that sub-exponential forms () are preferred, with the best-fit value of 0.590.02 for the young pulsars and 0.70.15 for the millisecond pulsars. The broken power-law fit to the Crab pulsar data from Aliu et al. (2011) is plotted for scale. Note that only statistical uncertainties on the SED data points were used during the fitting and thus the uncertainty on the best-fit parameter values are likely underestimated. This figure is taken from McCann (2014).

Iv Discussion and Conclusion

Following a stacked analysis of 115 gamma-ray pulsars, with an average exposure of 4.2 yr per pulsar, we find no evidence of cumulative emission above 50 GeV. Stacked searches exclusive to the young pulsars, the millisecond pulsars, and several other promising sub-samples also return no significant excesses above 50 GeV. Any average emission present in the entire pulsar sample is limited to be below 7% of the Crab pulsar in the 56-100 GeV band. The average flux limits presented in Table 1 are roughly 3 times lower than the best flux limits achieved in dedicated individual pulsar analyses done in the 2PC in the 56-100 GeV band.

One should note that a limit on the average flux from 115 pulsars at 7% of the Crab pulsar level is consistent with, for example, a scenario in which all 115 pulsars emit at 7% of the Crab pulsar level. It is also consistent with a scenario in which 8 pulsars emit at 100% the level of the Crab pulsar and the remaining 107 pulsars have zero emission. Therefore this analysis does not exclude the possibility of finding several pulsars which are as bright as the Crab pulsar above 50 GeV, or several dozen which are ten times dimmer. It does, however, constrain the average flux from the ensemble, and therefore for every individual pulsar detected above this flux limit, the average emission from the remaining pulsars is constrained to be further below the limit.

In the 100 MeV to 50 GeV energy range we find that the average SEDs returned from the young pulsar and millisecond pulsar stacking analyses are very similar in shape and are generally compatible with a power law times a sub-exponential cutoff. Abdo et al. (2010) and Celik & Johnson (2011) have shown that a sub-exponential cutoff function approximates a superposition of exponential cutoffs, thus the appearance of a sub-exponential cutoff in the ensemble SED is to be expected within a curvature radiation model. We note, however, that the highest energy spectral point is higher than the best fit sub-exponential cutoff function at the 2.4 level in both the young pulsar and millisecond pulsar cases. This cannot be taken as strong evidence for a non-exponentially-suppressed pulsar emission component aggregating in the stacked analysis, however, the available data cannot rule it out beyond the level of the limits shown in Figure 3 and Table 1.

Beyond this work, improvements can be made using the forthcoming Fermi-LAT pass-8 data release which will improve the Fermi-LAT acceptance by 25% at 100 GeV (Atwood et al., 2013). Improvements to this stacking analysis can also be made by employing a likelihood framework to stack the sources (see Ackermann et al. 2011 for example), rather than the simple On minus Off procedure described here. The flux sensitivity of any stacking analysis will, however, ultimately be bounded by the exposure of the Fermi-LAT. A future stacking analysis which doubles both the number of pulsars and the duration of observation used will increase the exposure by a factor of 4, indicating that future stacking analyses which do not yield detections may improve on the limits presented here by perhaps one or two orders of magnitude.

A more detailed account of the stacking analysis methods and results of this study can be found in McCann (2014).

This research is supported in part by the Kavli Institute for Cosmological Physics at the University of Chicago through grant NSF PHY-1125897 and an endowment from the Kavli Foundation and its founder Fred Kavli. I am grateful to the VERITAS collaboration for their support. In particular, I am very grateful to Pat Moriarty, David Hanna, Nepomuk Otte, Benjamin Zitzer, Jeremy Perkins, Scott Wakely, Nahee Park, Jeffery Grube and Christopher Bochenek for very helpful discussions on aspects of the analysis presented here.


  1. At this symposium the MAGIC collaboration presented evidence indicating that the power-law spectrum of the Crab pulsar may extend to TeV energies. See http://fermi.gsfc.nasa.gov/science/mtgs/symposia/2014/abstracts/185
  2. In more realistic models the acceleration efficiency is expected to be a few percent to a few tens of percent (Lyutikov et al., 2012).
  3. http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/aperture_photometry.html
  4. https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars
  5. A total of 117 pulsars are listed in the 2PC, however, the Crab pulsar was excluded from this analysis since we are investigating whether high-energy Crab-pulsar-like emission is seen in other pulsars. Further, PSRJ2215+5135 was also excluded from the study since no Off phase region was listed for this source in the 2PC.
  6. The amount of data analyzed here depends entirely on the availability and validity of pulsar spin-down timing solutions used for phase-folding.
  7. The Crab pulsar was excluded from all of these sub-sample stacking analyses. The parameter values listed in the 2PC catalog were used in all cases.


  1. Abazajian, K. N., & Kaplinghat, M. 2012, Phys. Rev. D , 86, 083511
  2. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ , 713, 154
  3. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010b, ApJS , 187, 460
  4. Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS , 208, 17
  5. Ackermann, M., Ajello, M., Albert, A., et al. 2011, Physical Review Letters, 107, 241302
  6. Ackermann, M., Ajello, M., Allafort, A., et al. 2013, ApJ , 765, 54
  7. Aharonian, F., Akhperjanian, A. G., Barrio, J. A., et al. 1999, A&A , 346, 913
  8. Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2004, A&A , 421, 529
  9. Aharonian, F. A., Bogovalov, S. V., & Khangulyan, D. 2012, Nature , 482, 507
  10. Akerlof, C. W., Breslin, A. C., Cawley, M. F., et al. 1993, A&A , 274, L17
  11. Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJ , 674, 1037
  12. Aliu, E., The VERITAS Collaboration et al. 2011, Science, 334, 69
  13. Aliu, E., Archambault, S., Archer, A., et al. 2015, ApJ , 800, 61
  14. Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514
  15. Bednarek, W. 2012, MNRAS, 424, 2079
  16. Celik, O., & Johnson, T. J. 2011, American Institute of Physics Conference Series, 1357, 225 utoffs
  17. Du, Y. J., Qiao, G. J., & Wang, W. 2012, ApJ , 748, 84
  18. Harding, A. K., Stern, J. V., Dyks, J., & Frackowiak, M. 2008, ApJ , 680, 1378
  19. Helene, O. 1983, Nuclear Instruments and Methods in Physics Research, 212, 319
  20. Hirotani, K. 2001, ApJ , 549, 495
  21. Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
  22. Leung, G. C. K., Takata, J., Ng, C. W., et al. 2014, ApJ , 797, LL13
  23. Li, T.-P., & Ma, Y.-Q. 1983, ApJ , 272, 317
  24. Lyutikov, M., Otte, N., & McCann, A. 2012, ApJ , 754, 33
  25. Lyutikov, M. 2012, ApJ , 757, 88
  26. McCann, A. 2014, arXiv:1412.2422
  27. Neshpor, Y. I., Stepanyan, A. A., Zyskin, Y. L., et al. 2001, Astronomy Letters, 27, 228
  28. Pétri, J. 2012, MNRAS, 424, 2023
  29. Singh, B. B., Chitnis, V. R., Bose, D., et al. 2009, Astroparticle Physics, 32, 120
  30. Vishwanath, P. R., Sathyanarayana, G. P., Ramanamurthy, P. V., & Bhat, P. N. 1993, A&A , 267, L5
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description