---
abstract: 'We perform spectral analyses of four bright $\gamma$-ray pulsars: PSR J0007+7303, Vela, Geminga and PSR J2021+3651 concentrating on the high-energy tails, defined as emission above 10GeV. The two competing models of pulsar $\gamma$-ray emission predict qualitatively different spectra well above the break energy: curvature emission predicts an exponential cut-off in the spectra, while the inverse-Compton scattering mechanism favors a power-law. We perform fits to the phase-averaged spectral energy distributions for each of the four pulsars. We find that in all cases both models fit the data equally well—the present data set does not allow any firm claim to be made about the shape of the spectra above 10GeV. In no case is the power-law fit or exponential cut-off fit significantly preferred over the other. The Crab pulsar remains the only known pulsar in which the power-law fit is clearly preferred over the exponential cut-off.'
author:
- 'G. T. Richards, M. Lyutikov,'
bibliography:
- '/Users/maxim/Home/Research/BibTex.bib'
title: 'High energy tails of pulsar gamma-ray emission'
---

Introduction
============

At present more than 200 pulsars[^1] are known to emit $\gamma$-rays in the high-energy (HE; E $\gtrsim$ 100MeV) band, and this number is steadily growing [@2013ApJS..208...17A]. Understanding the emission mechanism of $\gamma$-ray pulsars and its relationship to radio pulsar emission is one of the most pressing questions in high-energy astrophysics.

Since the launch of the [[*Fermi*]{}]{} satellite, the $\gamma$-ray data unambiguously point to an outer-magnetospheric origin of the high-energy emission. The observed spectra lack the super-exponential cut-off expected in “polar cap” models [@2008Sci...322.1221A], and there are far more radio-quiet $\gamma$-ray pulsars than are predicted in the polar cap scenario. There are a few independent [*geometrical*]{} models of pulsar high-energy emission that assume enhanced particle acceleration and radiation along the last open field lines of a magnetic dipole [@1995ApJ...438..314R; @2010ApJ...715.1282B]. Though [*geometrical*]{} models of the $\gamma$-ray emission from pulsars have been quite successful, the underlying physics of particle acceleration and emission remains an open question. In modeling the electrodynamics of the pulsar gap, it is typically assumed that the primary emission mechanism is curvature radiation [[*e.g., *]{} @1996ApJ...470..469R; @Baring99]. Our current theoretical models of pulsar radiation are not yet able to produce solid predictions. Indeed, there are several unanswered questions: the global kinetic structure of the magnetosphere, the location of the gaps, the current closure, and the location of pair production [@2014ApJ...795L..22C; @2015MNRAS.448..606C; @2015MNRAS.449.2759B]. Also, the process is non-local—photons emitted at one place can produce an avalanche on different field lines, so the gap may not self-close even though it produces a cascade.

The VERITAS Collaboration detected pulsed emission from the Crab pulsar at energies above 100GeV [@2011Sci...334...69V]. This discovery was confirmed by MAGIC, who have recently claimed the detection of the pulsed emission up to 1.5TeV . Most importantly, the photon spectrum of the pulsed emission between 100MeV and 400GeV is best described by a broken power law and is statistically preferred over a power law with an exponential cut-off.

Conventionally, $\gamma$-ray emission from pulsars is attributed to incoherent curvature radiation. The detection of the Crab pulsar by VERITAS [@2011Sci...334...69V] implies the importance of inverse-Compton (IC) scattering for the production of $\gamma$-ray photons. This follows from a very general argument for an upper limit on the possible energies of curvature photons [following a similar approach applicable to synchrotron emission; @1996ApJ...457..253D; @2010MNRAS.405.1809L]. Assuming that the radius of curvature $R_c $ of the [[magnetic field]{}]{} lines is a fraction $\xi$ of the light cylinder radius $R_c = \xi R_{LC}$ and balancing radiative losses with acceleration by the [[electric field]{}]{} (which is a fraction of the [[magnetic field]{}]{}, $E=\eta B$), one can find [*the maximum energy of curvature emission within the Crab pulsar [magnetosphere]{}*]{}: $
\epsilon_{br} 
 =
(150\, {\rm GeV}) \, \eta^{3/4} \sqrt{\xi} .
$ [@2012ApJ...754...33L]. The possibility that the emission above the break energy is produced as a tail to the curvature emission is excluded by the fact that the spectrum of the Crab pulsar reported by MAGIC and VERITAS [*does not show the exponential cut-off*]{} indicative of radiation reaction-limited curvature emission. Phase-dependent variations of the break energy at few GeVs should not affect the shape of the spectrum at sufficiently high energies, [*e.g., *]{}above 10GeV.

Importantly, it is the high-energy tails that carry information about the emission mechanism and thus deserve special consideration. The two competing models of pulsar $\gamma$-ray emission predict qualitatively different spectra well above the break energy (curvature emission predicts an exponential cut-off, while inverse-Compton scattering predicts a power-law).

One of the problems of the broadband double power-law fits, commonly used for a typical [*Fermi*]{}-LAT spectrum of a bright $\gamma$-ray pulsar, is that most of the errors accumulate due to the [*arbitrary*]{} parametrization of the spectral near the peak, where the two power laws connect. A different parametrization of this spectral part can produce substantially different fits [@2012ApJ...757...88L].

[@2013MNRAS.431.2580L] showed that the broadband spectrum of the Crab pulsar from UV to very high-energy $\gamma$-rays—nearly ten decades in energy—can be reproduced within the framework of the cyclotron self-Compton model [see also @2015ApJ...811...63H]. Emission is produced by two counter-streaming beams within the outer gaps, at distances above $\sim$20 NS radii. The outward moving beam produces UV–X-ray photons via Doppler-boosted cyclotron emission, and GeV photons by inverse-Compton scattering the cyclotron photons produced by the inward going beam. The scattering occurs deep in the Klein-Nishina regime, whereby the IC component provides a direct measurement of the particle distribution within the magnetosphere. The required plasma multiplicity is high, on the order of $10^6$ to $10^7$, but is consistent with the average particle flux injected into the pulsar wind nebula. Klein-Nishina reduction in the scattering cross-section (and the corresponding reduction of the electron energy loss rate) allows the primary leptons to be accelerated to very high energies with hard spectra. The secondary plasma is less energetic but more dense and has approximately the same energy content as the primary beam, producing cyclotron self-Compton radiation. The cyclotron component has a broad peak in the UV–X-ray range, while the IC component extends to hundreds of GeV. Below $\sim$1GeV the curvature emission from the primary beam can contribute a substantial flux fraction. The IC emission from the primary beam extends well into the TeV regime but has proven difficult to detect by Cherenkov telescopes.

How special is the Crab in terms of the high-energy power-law tail? [@2012ApJ...757...88L] analyzed the high energy emission from Geminga (using post-processed data) and concluded that the high energy part of the spectrum is better fit with a power-law shape than with an exponential cut-off. The non-detection by VERITAS [@2015ApJ...800...61A] is consistent with the fairly steep power-law slope inferred by [@2012ApJ...757...88L]. [@2014ApJ...797L..13L] have reported the detection of pulsed emission from Vela $> 50$GeV, noting that the observed spectrum above 10GeV of the Vela pulsar is harder than a simple exponential function predicted in the exponential cut-off spectral shape of the tail. [@2015ApJ...804...86M] has performed a statistical analysis of high energy tails using a stacked analysis of 115 [*Fermi*]{}-LAT-detected pulsars, finding no significant stacked excess at energies above 50GeV.

In this paper we present the results of a complementary analysis of high-energy pulsar tails, defined as emission above 10GeV. We chose four bright $\gamma$-ray pulsars: PSR J0007+7303, PSR J0633+1746 (Geminga), PSR J0835–4510 (Vela) and PSR J2021+3651, and we perform detailed spectral fits for each pulsar.

Fermi-LAT Data Analysis & Results
=================================

The [*Fermi*]{}-LAT is an electron-positron pair-conversion telescope sensitive to gamma-ray photons with energies between 20MeV and 300GeV. It has a FoV of $\sim2.5$sr and attains full-sky coverage approximately every three hours. For a complete description of the instrument, see [@2009ApJ...697.1071A; @2012ApJS..203....4A].

To obtain the spectral parameters for each of the four pulsars, the [*Fermi*]{}-LAT `Science Tools v10r0p5` with Pass 8 reprocessed instrument response functions and the standard quality cuts described in [@2012ApJS..199...31N] are used. Roughly 7yr of “source”-class events with energies between 100MeV and 100GeV collected between 4 August 2008 and 12 August 2015 within a $20^\circ$ region-of-interest (ROI) of the 2PC location of each pulsar are processed with the maximum likelihood fitting routine.

The spectral analysis for each pulsar except Vela presented here follows the same steps outlined in the second [*Fermi*]{}-LAT pulsar catalog [2PC @2013ApJS..208...17A]. To generate spectral energy distributions of the [*Fermi*]{}-LAT data for each pulsar, a binned maximum likelihood analysis is done in logarithmically spaced energy bins spanning the range 100MeV to 100GeV. Spectral models for all sources in the 3FGL catalog [@2015ApJS..218...23A] in the ROI in addition to the galactic and isotropic diffuse backgrounds (gll\_iem\_v06.fits, iso\_P8R2\_SOURCE\_V6\_v06.txt) are included in the likelihood fitting. Sources within a circle $4^{\circ}$ in radius are modeled as power laws of the form $${dN \over dE} = N_0 {\left(E \over E_0\right)}^{\gamma},$$ where $N_0$ is the flux normalization, $E_0$ is fixed to 300GeV, and $\gamma$ is the spectral index. These sources have only their normalization parameter left free in the fitting routine. All other parameters are fixed except for the normalization on the galactic and isotropic diffuse background models. In each bin, the pulsar is modeled as a point source with a simple power-law spectral shape with a fixed spectral index of $-2$. An iterative spectral binning is performed starting at three bins per decade in energy and increasing to four bins per decade, five, and so on. The spectrum containing the greatest number of bins above 10GeV is kept for subsequent fitting, while the others are discarded.

The [*Fermi*]{}-LAT data analysis for the Vela pulsar is done differently than described above due to three nearby bright diffuse sources (Vela X, Vela Jr., and Puppis A) that must be taken into account in order to achieve the greatest sensitivity. To extract spectral parameters for the diffuse sources, the data are first phase-folded with the the `Tempo2` pulsar timing package [@2006MNRAS.369..655H] using a publicly available timing solution from a webpage maintained by M. Kerr[^2]. Phase-folding the data allows selection of the off-pulse region to remove the contribution from the Vela pulsar, which allows more robust likelihood fitting for the nearby diffuse sources. Since the timing solution is not valid for the full 7yr data set, only data recorded by the [*Fermi*]{}-LAT during the period of validity of the timing solution are used in the following analysis. First, the off-pulse phase region between 0.8 – 1.0 is selected to remove the Vela pulsar. Subsequently, all normalization parameters in the model are scaled to 20% of their nominal values, and the Vela pulsar model is removed. All sources within $4^{\circ}$ of the Vela pulsar 2PC location are modeled as power laws with free normalization parameters, including the three diffuse sources and the galactic and isotropic diffuse components. Likelihood analyses are done in logarithmically spaced energy bins to allow computation of the normalization parameters for each source in each bin. The normalization parameters are then scaled by a factor of 5, and the spectral model for the Vela pulsar is added. The normalization parameters for the three diffuse sources are then fixed to the previously computed values in each energy bin, and the Vela pulsar SED is generated in the same way as described in the previous paragraph.

Phase-averaged spectral energy distributions for each pulsar are shown in Figure \[seds\]. Each SED is fitted with both a power-law and an exponential cut-off above 10GeV to test the preferred spectral shape. The choice of 10GeV is an arbitrary starting point for probing the highest energies of the spectra. However, 10GeV is above a few times the break energy for each pulsar reported in the 2PC, which helps to reduce contamination by the harder portion of the spectra at lower energies. Spectral indices for the power-law fits and fit probabilities are summarized in Table \[tab:fit\_summary\].

[1.0]{}[CCCCC]{} Pulsar & Power-law spectral index above 10GeV & Power-law fit chi-square / ndf (probability) & Exponential cut-off fit chi-square / ndf (probability) & Notes\
PSR J0007+7303 & $-3.69 \pm 0.016$ & 5.26 / 3 (0.15) & 0.95 / 3 (0.81) &\
PSR J0633+1746 & $-5.12 \pm 0.16$ & 5.54 / 2 (0.06) & 4.87 / 2 (0.09) & Geminga\
PSR J0835–4510 & $-4.54 \pm 0.08$ & 14.1 / 5 (0.02) & 8.06 / 5 (0.15) & Vela\
PSR J2021+3651 & $-4.73 \pm 0.54$ & 1.10 / 1 (0.29) & 2.40 / 1 (0.12) &\

![image](combined_4panel_SED_v4.pdf){width="6.6in"}

The fit probabilities for either the power-law fit or the fit with an exponential cut-off are comparable. Therefore, one has to conclude that both models describe the data equally well and, therefore, both models are equally preferred. The probability of a fit to the Vela data with a power law is the worst of all four power law fits. In the case of Vela, the data points at the lowest energies included in the fit contribute the most to the $\chi^2$. In that energy range the spectrum shows clear curvature, highlighting the problem of arbitrarily choosing 10GeV as the energy threshold for the fit. It also emphasizes that the sensitivity of the method strongly depends on how far away from the break in the spectrum the fit begins.

In order to alleviate the subjective choice of where to start the fit, the spectra are also fitted with a power law plus sub-exponential cutoff (PLSEC) of the form $${dN \over dE} = N_0 {\left(E \over E_0\right)}^{\gamma_1} \textrm{exp}\left(-{\left(E \over E_c\right)}^{\gamma_2}\right),$$ where $E_0$ is fixed to 300GeV, and $N_0$, $E_c$, $\gamma_1 $ and $\gamma_2$ are free parameters, and a smooth, broken power law (SBPL) of the form $${dN \over dE} = N_0 {\left(E \over E_0\right)}^{\gamma_1}\left( 1 + {\left(E \over E_b \right)}^{\gamma_1 - \gamma_2 \over \beta} \right)^{-\beta},$$ where $E_0$ is again fixed to 300GeV, and $N_0$, $\gamma_1$, $\gamma_2$, $E_b$, and $\beta$ are all left free in the fit. The parameters of the fit functions are constrained by fitting them to data points below 10GeV. For both the PLSEC and the SBPL, including the minimum number of data points below 10GeV to constrain the parameters of the fit yields acceptable fit probabilities. However, including more data points at lower energies in the fits, the fits fail with probabilities quickly sinking below $10^{-5}$. For all energy ranges for the fits, the PLSEC clearly under-predicts the flux above 10GeV, and the SBPL law clearly over-predicts the data—neither of the fitted functions match the data above 10GeV.

We note that the upper limit at 48.4GeV for PSR J0007+7303 is a factor of $\sim$3.5 below the power-law fit, but it still lies within 1.2$\sigma$ of the fit taking the error on the fit into account. Therefore, the upper limit at 48.4GeV is not inconsistent with the power-law scenario.

Discussion
==========

We have performed an analysis of the high energy tails, defined as emission above 10GeV, of four of the brightest [*Fermi*]{}-LAT-detected $\gamma$-ray pulsars. Inn our approach, the binning was adjusted iteratively starting at three bins per decade and increasing, keeping the spectrum with the most bins above 10GeV for subsequent fitting. In the second [*Fermi*]{}-LAT pulsar catalog [2PC @2013ApJS..208...17A] the binning was arbitrary, mostly tailored to weaker sources. As a result, for the bright sources considered in this work, the binning was not optimal.

The fit results for all four of the pulsars analyzed in this study do not allow any firm claim to be made about the shape of the spectra above 10GeV. In no case is the power-law fit or exponential cut-off fit significantly preferred over the other, though the exponential cut-off is marginally preferred in the case of PSR J0007+7303 (which is the dimmest pulsar of all considered in this work). Previously, in the case of Geminga, [@2012ApJ...757...88L] used post-processed data, with the standard binning, while in the present work binning was adjusted to produce maximum likelihood of a fit. One of the unquantified biases in our analysis is the choice of the lower energy cut-off, taken here to be 10GeV. This choice is motivated by our attempt to avoid fitting the roll-off region near the break energies.

Alternatively we chose the lower cut at three times the peak energy. In the sample, PSR J0007+7303 has the highest peak energy of $4.7$ GeV and is the weakest of all, so a choice of higher energy cut reduces the statistical significance. The break energies in Vela and PSR J2021+3651 are both $3.0$ GeV. The only exception is Geminga, with the break energy of $2.2$ GeV. Our work leaves open the question of the shape of the high-energy gamma-ray spectral tails of these four pulsars.

[^1]: 

[^2]: ``
