{'fraction_non_alphanumeric': 0.03693501112898605, 'fraction_numerical': 0.03323269497763184, 'mean_word_length': 4.213465073529412, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PSR J0007+7303 in the CTA1 SNR: New Gamma-ray Results from Two Years of Fermi-LAT Observations

One of the main results of the Fermi Gamma-Ray Space Telescope is the discovery of {\gamma}-ray selected pulsars. The high magnetic field pulsar, PSR J0007+7303 in CTA1, was the first ever to be discovered through its {\gamma}-ray pulsations. Based on analysis of 2 years of LAT survey data, we report on the discovery of {\gamma}-ray emission in the off-pulse phase interval at the ~ 6{\sigma} level. The flux from this emission in the energy range E \geq 100 MeV is F_100 = (1.73\pm0.40)\times10^(-8) photons/cm^2/s and is best fitted by a power law with a photon index of {\Gamma} = 2.54\pm0.14. The pulsed {\gamma}-ray flux in the same energy range is F_100 = (3.95\pm0.07)\times10^(-7) photons/cm^2/s and is best fitted by an exponentially-cutoff power-law spectrum with a photon index of {\Gamma} = 1.41 \pm 0.23 and a cutoff energy E_c = 4.04 \pm 0.20 GeV. We find no flux variability neither at the 2009 May glitch nor in the long term behavior. We model the {\gamma}-ray light curve with two high-altitude emission models, the outer gap and slot gap, and find that the model that best fits the data depends strongly on the assumed origin of the off-pulse emission. Both models favor a large angle between the magnetic axis and observer line of sight, consistent with the nondetection of radio emission being a geometrical effect. Finally we discuss how the LAT results bear on the understanding of the cooling of this neutron star.


Introduction
Theşupernova remnant (SNR) (G119.5+10.2) is a composite SNR characterized by a large radio shell enclosing a smaller pulsar wind nebula (PWN). Pineault et al. (1993) derived a kinematic distance of 1.4 ± 0.3 kpc based on associating an H I shell found to the northwestern part of the remnant to the remnant itself. Observations in X-rays with ASCA and ROSAT revealed a central filled SNR with emission extending to the radio shell in the south and southeast as well as in the north and northwest radio-quiet regions (Seward et al. 1995). X-ray observations with ROSAT also revealed the point source RX J0007.0+7302 (Seward et al. 1995). X-ray observations with Chandra revealed a compact PWN and a jetlike structure (Halpern et al. 2004). Radio and X-ray characteristics of imply an age in the range 5000-15000 years (Pineault et al. 1993;Slane et al. 1997Slane et al. , 2004. Slane et al. (2004) estimated for the age of the SNR a value of 1.3 × 10 4 d 1.4 yr (where d 1.4 is the distance in units of 1.4 kpc) which is in good agreement with the spin-down age estimate of 14,000 yr from γ-rays (Abdo et al. 2008). The offset of the point sourcefrom the geometrical center of the radio SNR allows for the estimate of the transverse velocity of the point source which Slane et al. (2004) estimated to be ∼450 km s −1 .
Prior to the launch of Fermi the EGRET γ-ray source 3EG J0010+7309, which lies within the boundaries of the radio SNR, showed the characteristics of a pulsar (Brazier et al. 1998;Halpern et al. 2004). Mattox et al. (1996b) discussed this source as a potential candidate for a radio-quiet γ-ray pulsar. A search for γ-ray pulsations using EGRET data didn't reveal the pulsar (Ziegler et al. 2008). The characteristics of this source in X-rays also pointed to it as a pulsar (Halpern et al. 2004). In fact Halpern et al. (2004) called the source the "pulsar" although no pulsations were detected from this source using FFT searches on ∼ 26 ks of XMM-Newton data (Slane et al. 2004) or in radio with the Green Bank Telescope (Halpern et al. 2004). Very deep searches for a counterpart forin optical and radio resulted only in upper limits (Halpern et al. 2004). By correlating X-ray images from Chandra with those in the optical waveband of theŗegion, Halpern et al. (2004) gave the most accurate position of this source to date ((J2000.0) 00 h 07 m 1. s 56,+73 • 03 ′ 08.1 ′′ ) with an accuracy of ∼ 0 ′′ .1. The discovery of the pulsar didn't come until the launch of Fermi. During its commissioning stage the LAT discovered a 315.87 ms pulsar at the location of RX J0007+7303 (Abdo et al. 2008). The discovery of the pulsar in γ-rays prompted a long, ∼ 130 ks, XMM-Newton observation to search for pulsations from this source (PI: Caraveo 2008, ObsID: 06049401). Using this XMM-Newton data set along with the timing model from the LAT (Abdo et al. 2009b) X-ray pulsations from this source in X-rays were finally detected (Lin et al. 2010;Caraveo et al. 2010).
This paper reports further analysis on the LAT data and related observations of the ŞNR and its γ-ray pulsar, extending the initial report of the pulsar discovery (Abdo et al. 2008). The new developments contribute to a characterization both of the pulsar and its relation to the associated extended source, exploiting two years of data now accumulated by Fermi. New developments fall in three distinct areas. First, timing analysis of the pulsar spanning the two years shows a glitch near the middle of the interval at MJD 54952.652 (1 May 2009), with relative change in pulse frequency ∆ν/ν ∼ 6 × 10 −7 . This permits careful comparison of the pulsar characteristics before and after the glitch. Second, we present new γ-ray results concerning the relationship of the pulsar to the surrounding remnant, the first detection of an extended source in the off-pulse emission. Finally, we summarize the multi-wavelength picture of the source. From this perspective the outstanding characteristic of the neutron star in CTA1 is that it appears cool for its inferred age. We reconsider both the temperature and the age estimates, and then relate this to pulsar characteristics established from timing LAT γ-rays. In earlier work (Halpern et al. 2004) the pulsar in ḩas been assigned an inferred mass exceeding 1.42 M ⊙ even though it is not in a binary, because the larger mass should accelerate cooling. Even from initial timing solutions it was clear that this pulsar has a very strong magnetic field, though weaker than that of a magnetar, hence it is atypical in two respects.

Gamma-Ray Observations and Data Analysis
We used 2 years of survey data collected with the LAT to study this source in γ-rays. The data set starts 2008 August 4 and ends 2010 August 4 (54682.68 -55412.65 MJD). This large data sample with high statistics, compared to the six weeks of observations used for the discovery paper (Abdo et al. 2008) and the six months used for the Fermi pulsar catalog paper (Abdo et al. 2010b) allows us to perform several key timing and spectral analyses not feasible in prior studies. In particular we study the phase-resolved spectra and the flux variability especially around the 2009 May 1 glitch, we search for off-pulse emission and build a precise timing model. Throughout this paper we used "diffuse" class photons events with the P6 V11 instrument response functions (IRFs) 1 (Atwood et al. 2009). To reject atmospheric γ-rays from the Earth's limb, we selected events with zenith angle < 100 • .

Timing Analysis
For the pulse timing analysis, we selected events with energies > 170 MeV that were reconstructed within 1.3 • from the Chandra location for RX J0007.0+7303 source (Halpern et al. 2004). These radius and energy cuts were selected to maximize the pulsed significance. Following the procedure described by Ray et al. (2011), we measured a total of 72 pulse times of arrival (TOAs), each with a integration time of about 10 days, referenced to the geocenter. Each TOA was determined using the unbinned maximum likelihood technique from Ray et al. (2011) using a template profile consisting of 2 gaussian components. The TOAs were then fit to a timing model using Tempo2 (Hobbs et al. 2006). The model, shown in Table 1, includes position, frequency (ν) and the first three frequency derivatives (ν,ν, ... ν ). We fitted for position since this will give the smallest timing residuals and will allow for the comparison of our timing position to that from Chandra. We find our timing position, shown in Table 1, to be consistent with the Chandra position of compact source RX J0007.0+7303 (00:07:01.56, 73:03:08.3; see Halpern et al. (2004)). In addition, a glitch on 2009 May 1 (MJD 54952.652) was included with an instantaneous step in ν andν at the glitch. The glitch epoch was chosen to produce a zero phase jump at the glitch. For this glitch we measure ∆ν/ν = 5.54(1) × 10 −7 . With a shorter dataset, Ray et al. (2011) could not be certain of the step inν at the glitch. Our extended observations allow us to be confident of the ∆ν/ν = 9.7(6) × 10 −4 at the glitch.
The timing model includes second and third frequency derivatives, which are required to obtain white residuals. This is presumably an indication of timing noise in this young pulsar. We note that for a braking index of 3, aν of 1.24 × 10 −23 s −3 is expected from the secular spin down of the pulsar. This accounts for about a third of the total measuredν. If one interprets the measuredν as being entirely due to the secular spin down of the pulsar, one obtains a braking index (n) of n = νν/ν 2 = 9.95.

Detection of Off-Pulse Emission
To search for any γ-ray emission present in the off-pulse part of the phase, we had first to determine accurately the definition of the off-pulse phase window. One might simply determine the off-pulse interval by eye. In that method the off-pulse starts when any apparent pulsed emission decreases to the levels of the background and ends when the pulsed emission resumes and an increase above background is seen. A less arbitrary method is to perform a likelihood analysis in small phase bins, gradually increasing the width of consecutive bins. For all of these width-varying bins only the left and right limits are changing while the center of the bin is fixed at what one believes to be the center of the off-pulse interval. A spectral shape is assumed for any off-pulse emission and the increase in signal as a function of phase bin width (∆φ) is analysed. We performed this analysis with gtlike where we selected a region of interest of 20 • and a source region of 30 • (see §2.3 for more details). In this analysis all phase bins were centered at φ = 0.89. In the case of absence of signal in off-pulse, one would expect the distribution of the test statistic (TS) values (Mattox et al. 1996a) to be centered around zero with no correlation with the width of the phase bins ∆φ. However, in the case of presence of off-pulse emission one would expect the significance to increase linearly with the increase in the width of the phase bins until one starts integrating photons from the pulsar itself where a sharp increase in significance is then seen. The off-pulse width is then defined as the point at which the sharp increase in the TS occurs. This is shown in Figure 1.
From the plot, one sees two ranges of data points, with a clear break in slope between the earlier range (first seven points) and the remainder. We fitted first order polynomials to each range and detected a clear significant change in slope, by about a factor of four. This is taken to be the onset of contamination from the pulsed signal. We therefore define the off-pulse interval to be φ ∈ [0.71 − 1.07] with a width of ∆φ = 0.36.
As can be seen from Figure 1 there is a significant detection of γ-ray emission in the off-pulse phase. At a TS of ∼ 40 this is the first detection of γ-ray emission in the off-pulse of PSR J0007+7303. An earlier Fermi -LAT survey of PWNe by Ackermann et al. (2011), using 16 months of LAT data, noted a candidate for an off-pulse emission from PSR J0007+7303 but it was below the detection threshold even though a wider phase window was used for the off-pulse. The present detection, using 2 years of LAT data but with a conservative phase window is unambiguous. Figure 2 shows a TS map of the off-pulse part of PSR J0007+7303. On the same Figure  we show ROSAT X-ray contours in black (Seward et al. 1995) and radio contours in green (Pineault et al. 1997). From the Figure one can see that 1) the γ-ray signal detected with the LAT is better correlated with the ROSAT X-ray emission than with the radio SNR and 2) that the off-pulse γ-ray emission is clearly extended.

Off-Pulse Extension Analysis
To check for a possible extension in the off-pulse γ-ray emission, we performed the likelihood analysis similar to that in §2.3.1 with the further addition of an extended disk template to describe the possible extension of the emission. Different angular sizes in the range 0.1 • -1.0 • and centroid position templates have been fitted. In each case the significance was compared to a simple point-like hypothesis. A disk of radius 0.6 • ± 0.3 • is favored, and a point-like hypothesis is excluded at the 95% confidence level. A template in the shape of an ellipse was also fitted to the off-pulse emission. We found no compelling statistical evidence in favor of this fit compared to the disk.

Spectral Analysis
Spectral analyses for this source were performed using the Fermi LAT maximumlikelihood Science Tool gtlike in its binned mode 2 . Fits were performed on a 14 • × 14 • region of the sky centered at the pulsar position selecting photons in the energy range 0.1 to 300 GeV. We used a model that included diffuse emission components as well as nearby γ-ray sources from the First Fermi -LAT γ-ray catalog (1FGL) (Abdo et al. 2010a) that fell within 19 • from the position of PSR J0007+7303 . The Galactic diffuse emission was modeled using the gll iem v02 P6 V11 DIFFUSE model and the isotropic background using the isotropic iem v02 P6 V11 DIFFUSE model 3 .
In performing the fit we fixed all the parameters of the sources that fell between 14 • and 19 • from PSR J0007+7303 to their values in the 1FGL catalog, and left free the normalization factor of all the sources within 14 • of PSR J0007+7303. All the non-pulsar sources were modeled with a power law as reported in the 1FGL catalog, while the two pulsars in the region of interest, PSR J0205+6449 and PSR J2229+6114, were modeled by a power law with exponential cutoff according to the data reported in the Fermi -LAT pulsar catalog (Abdo et al. 2010c).
To obtain Fermi -LAT spectral points we divided our sample into logarithmically-spaced energy bins (4 bins per decade starting from 100 MeV) and then applied the maximum likelihood method in each bin. For each energy bin, all point sources, including PSR J0007+7303, were modeled by a power law with fixed photon index. From the fit results we then evaluated the integral flux in each energy bin. If in an energy bin the source significance is lower than 3σ we have evaluated the 95% integral flux upper limit in that bin. This method does not take into account energy dispersion or correlations among the energy bins. To obtain the points of the spectral energy distributions (SEDs) we multiplied the flux in each bin by the spectrally weighted mean bin energy.

Off-pulse Spectrum
To quantify the off-pulse γ-ray emission we have assumed a point-like source, modeled with a power law at the pulsar position. We considered only events in the off-pulse phase interval [0.71 -1.07]. The fitted power law spectrum is given by: where for this fit N = 1.69 ± 0.40 stat ± 0.18 sys × 10 −8 photons cm −2 s −1 , γ = 2.54 ± 0.14 stat ± 0.05 sys with E min = 100 MeV and E max = 100 GeV. The estimated integral flux above 100 MeV is F 100 = 1.73 ± 0.40 stat ± 0.18 sys × 10 −8 photons cm −2 s −1 and the integral energy flux above 100 MeV is G 100 = 7.83 ± 1.43 stat ± 0.56 sys × 10 −12 erg cm −2 s −1 . There is no compelling statistical case in favor of a cutoff in the spectrum, suggesting that the emission is not magnetospheric in origin. We choose to model the off-pulse emission with a power law because it has less parameters and gives smaller statistical errors compared to an exponential law with a cutoff.
Systematics are mainly based on uncertainties in the LAT effective area derived from the on-orbit estimations, and are of ≤ 5% near 1 GeV, 10% below 0.1 GeV and 20% above 10 GeV (Abdo et al. 2009a). We therefore propagate these uncertainties using modified effective areas bracketing the nominal ones (P6 V11 DIFFUSE).

On-pulse Spectrum
To account for the off-pulse emission we used the results of the off-pulse fit, properly rescaled to the on-pulse phase interval, as a starting point for the pulsed emission analysis. In the model we considered two sources in the same position, one described as a power law with the spectral parameters fixed at the values found with the off-pulse fit and one described by a power law with exponential cutoff (PLEC) in the form: We set b = 1 for which the best-fit parameters are N 0 = (9.08 ± 0.20 stat ± 0.54 sys )× 10 −11 cm −2 s −1 MeV −1 , Γ = (1.41 ± 0.23 stat ± 0.03 sys ) and E c = (4.04 ± 0.20 stat ± 0.67 sys ) GeV. The integral flux above 100 MeV is F 100 = (3.95 ± 0.07 stat ± 0.30 sys )× 10 −7 photons cm −2 s −1 and the corresponding energy flux is G 100 = (4.41 ± 0.06 stat ± 0.5 sys )× 10 −10 erg cm −2 s −1 . We also fitted the on-pulse phase-averaged spectrum to different spectral models, power law and broken power law. Both models can be excluded at the > 5σ level compared to the PLEC model used above. We have also fitted the PLEC spectrum while leaving free the exponential index b. The b value obtained was lower then 1 but the overall fit did not improve, thus we adopt the simpler PLEC model for which b=1. Figure 3 shows the results of the on-pulse phase-averaged spectrum.

On-pulse Phase-resolved Spectrum
To explore the on-pulse phase-resolved spectrum, we divided the pulse profile in variablewidth phase bins, each containing 500 photons above 100 MeV. These bins were defined by selecting only those events within an energy-dependent radius of θ < Max(Min(R max , θ 68 ), 0.35 • ) around PSR J0007+7303: the minimum value of 0.35 • was selected in order to keep all high-energy photons, while a maximum radius, R max = 2 • , was introduced to reduce the background contamination at low energies. This choice of binning provides a reasonable compromise between the number of photons needed to perform a spectral fit and the length of phase intervals. It should be short enough to sample fine details on the light curve, while remaining comfortably larger than the rms of the timing solution. A binned maximum likelihood spectral analysis, similar to the analysis performed in §2.3.2, was performed in each phase bin with the exception of fixing the spectral parameters of all the nearby γ-ray sources and of the two diffuse backgrounds to the values obtained in the phase averaged analysis, rescaled for the phase bin width. Using the likelihood ratio test we can reject the power law model at a significance level greater than 5σ in each phase interval. Such a model yields a robust fit with a logarithm of the likelihood ratio greater than 150 in each phase interval. Figure 4 shows the evolution of the spectral parameters across PSR J0007+7303's rotational phase. In particular, the energy cutoff trend provides a good estimate of the high energy emission variation as a function of the pulsar phase. Table 2 summarizes the results of the spectral fit in each phase bin. Variations of both the photon index and the cutoff energy as a function of rotational phase are apparent. The photon index seems to show a rough symmetry centered at a phase point half way between the two peaks with the index increasing (softening) outwards. This is similar to what has been observed for PSR J1709−4429 (Abdo et al. 2010e). The cutoff energy evolves quite differently as a function of the rotational phase. It increases from a minimum value of 1.5 GeV at φ = 0.2 until reaching the first peak where it stays at a constant value of ∼ 4 GeV until the second peak is reached, after which it seems to fluctuate between 3 and 6 GeV before finally decreasing to its minimum value of 1.5 GeV at φ = 0.6.

Light Curve
We investigated the pulsar light curve in different energy bands by selecting events within 1.6 • of the pulsar position. This value was selected to maximize the signal-to-background ratio over the full energy range (E > 100 MeV). The energy-resolved light curve is shown in Figure 5. The top panel in the Figure shows the folded light curve for energies above 100 MeV while the rest of the panels show the light curve in exclusive energy bands. The dashed horizontal line shown in the top panel shows the estimated level of the background due to diffuse emission. This background estimate of 194 ± 5 counts/bin was obtained by simulating two years of data. We used the LAT Science Tool gtobssim and used for the input model the best fitted model from §2.3.2 but with PSR J0007+7303 and the off-pulse component ( §2.3.1) removed from the input model. The light curve shows two distinct peaks. We fitted the light curve by a double Gaussian for which the first peak and second peak are located at φ = 0.303 ± 0.002 and φ = 0.484 ± 0.002 respectively. The separation between the means of the two peaks is 0.181 ± 0.003 in phase. As can be seen from the Figure there is a significant evolution in the counts ratio of the two peaks P1/P2.

Flux Variability Analysis
To check for long-term stability in the flux we performed likelihood analysis similar to that in §2.3 but in 8-day time bins. Figure 6 shows the resulting fluxes. The length of the time bin was selected to allow for the accumulation of enough statistics to guarantee a good likelihood fit. To look for flux variability from the source we adopt the method outlined in (Abdo et al. 2010a). The source shows no variability on this time scale. Similar analyses were performed for 16, 32, and 64-day time bins, no significant modulation was found.
To check for any change in the spectrum of the pulsar due to the glitch we split the data in two bins around the glitch. The pre-glitch epoch spans the time range 54682.68 -54952.652 (MJD), while the post-glitch data spans the time range 54952.652 -55412.65 (MJD). We performed a likelihood analysis similar to that in §2.5 in these two time bins. In Table 3 we show the spectral fits for the pulsar in these two epochs. The flux and spectral parameters are in good agreement for the two epochs. No change in integral flux above 100 MeV is seen.
Recent variability detected in the Crab appears to come from the Nebula, not the pulsar (Abdo et al. 2011). Since PSR J0007+7303 is among the younger pulsars and has a complex PWN/SNR associated with it, it is reasonable to ask whether analogous variability occurs in this source. While §2.2 has shown strong evidence for an off-pulse component, tentatively extended, there are not enough statistics to explore variability in this component at all, and certainly not on the time scales detected in the Crab.

Geometrical Constraints from Light Curve Modeling
The pulsar emission mechanism is not well understood, and the distributions of pulsar magnetic inclination angles, efficiencies, and other physical characteristics are largely unknown. There are several competing emission models, of which the resulting pulse profiles depend on the geometry of the system, defined by the emission zone location and size, observer viewing angle, and magnetic inclination angle. Fitting the profiles of these models to observed light curves can provide insight on the true emission and viewing geometries, for example leading to better constraints on luminosity and efficiency through calculation of the flux correction factor f Ω (equation 4 of Watters et al. (2009)). The geometry of a given system likely also contributes to the detectable presence or absence of radio and perhaps X-ray emission. For example, if the inclination and observer angles are very different, we could see a γ-ray-only pulsar, as the narrower radio beam may not cross our line of sight.
To determine the geometry of the CTA 1 pulsar and distinguish between emission models, we compared the LAT light curve of PSR J0007+7303 with the predicted light curves from geometrical representations of two standard high-energy pulsar emission models, the outer gap (OG) (Romani & Yadigaroglu 1995) and slot gap (SG, also referred to as the two-pole caustic or TPC) (Muslimov & Harding 2004) models. These models were considered within the context of the vacuum retarded dipole magnetic field.
The light curves were simulated as in Dyks et al. (2004), with the geometry modified to represent each emission region. The OG lies along the last open field lines between the null charge surface and the light cylinder, R lc = c/Ω, with r cyl ≤ 1.0R lc as large as possible for each geometrical configuration. The SG extends from the surface to r cyl = 0.95R lc . In both geometries, the maximum emission altitude r max was treated as a model parameter. We note that our representation of the SG differs from that of the TPC model in Romani & Watters (2010), as we allow emission out to much higher altitudes (their emission is cut off at r cyl = 0.75R lc ). The TPC geometry is therefore treated as a subset of the SG, rather than being considered as a completely separate model. The SG model has emission throughout the gap, while the OG emission occurs only along the field line at the gap's innermost edge. In the simulations, the vacuum retarded dipole field is assumed in the observer's frame and then transformed to the co-rotating frame (CF). Photons are emitted tangent to the field in the CF, prior to the aberration calculation. A constant emissivity is assumed along the field lines in the CF. The special relativistic aberration leads to a bunching in pulse phase, or caustics, on the trailing side of the pulse, producing peaks in the light curvyes. For a given inclination angle α, gap width w in open volume units (r ovc , as in Dyks et al. (2004)), and maximum emission radius r max in units of R lc , the code produces light curves at all observer angles ζ.
The LAT light curve of PSR J0007+7303 has 32 bins and a background count level of 195 or 246 counts/bin, depending on whether the off-peak emission discussed in §2.2 is (respectively) magnetospheric or due to a PWN. We compared this light curve with light curves simulated from a representative sample of all possible geometries and rebinned to match the observed light curve. Our models had five free parameters, α, ζ, w, r max , and a phase shift ∆φ introduced in order to best match the LAT light curve. We considered values of α between 0 • and 90 • and ζ between 0.5 • and 89.5 • , each with resolution of 1 • ; values of width 0 ≤ w ≤ 0.3 with resolution 0.01 r ovc ; and maximum emission radii 0.7 ≤ r max ≤ 2 for the SG and 0.9 ≤ r ≤ 2 for the OG, with resolution 0.1 R lc in both cases. For each model, the emission altitude is limited by Min(r max , r cyl )-r max therefore may be larger than r cyl , but emission ceases at the cylindrical radius if not at the maximum emission radius. The phase shift 0 • ≤ ∆φ ≤ 360 • is arbitrary because PSR J0007+7303 is radio-quiet within the flux limits achieved thus far by radio telescopes; were the pulsar radio-loud, the shift would be constrained to be at most equal to the phase lag between the radio and γ-ray peaks (Dyks et al. 2004).
To find the best-fit parameters of each model, we used a Markov Chain Monte Carlo (MCMC) maximum likelihood routine as described in Verde et al. (2003) that explored the parameter space. We calculated the χ 2 between the LAT light curve and model light curve and used Wilks' Theorem, ∆ ln(L) = −∆χ 2 /2, to guide the Markov chains toward regions of high likelihood. We also used the ∆χ 2 test to find 3σ confidence intervals on each parameter. For this pulsar, the best-fit light curve has a very large χ 2 . This is expected, as a) the pulsar is bright and its light curve has relatively small error bars and b) the light curve simulations result from geometrical models of possible but simplified pulsar magnetospheres, rather than from a well-understood physical model. When finding confidence intervals in steps of ∆χ 2 from a large initial χ 2 value, the intervals will appear artificially small. We therefore re-scale the χ 2 values found for all sets of parameters in the MCMC chains so that the minimum reduced χ 2 = N dof = 27 (the scale factor is therefore N dof /χ 2 min ). The confidence intervals are then found with these modified χ 2 values. In this way, we find parameter ranges that may give rise to the high-energy light curve of PSR J0007+7303.
For the case where the off-peak emission is assumed to be magnetospheric (corresponding to a background level of 195 counts/bin), the best-fit parameters are (α, ζ, w, r, ∆φ) = (6 • , 74.5 • , 0.04r ovc , 1.0 R lc , 2 • ) with χ 2 /27 = 22.5 for the OG model and (8 • , 69.5 • , 0 r ovc , 1.9 R lc , -4 • ) with χ 2 /27 = 11.3 for the SG model. Almost identical parameters are found for the case where the off-peak emission is assumed to originate from a PWN. The confidence intervals, f Ω values, and reduced χ 2 for the best fit in each interval are given in Table 4. The absolute best fit is found with the OG model using the higher background level. Figure 7 shows the LAT light curve superposed with the best-fit model light curves, as well as reduced χ 2 contours in α and ζ with w and r fixed at the best-fit values. To be complete, we also fit the light curve with the TPC model, the results of which are not shown here; we find this lower-altitude subset of the SG model cannot reproduce the sharp double peaks with the correct spacing anywhere in parameter space. We therefore find that both outer magnetosphere models produce light curves consistent with that which is observed. For magnetospheric off-peak emission, the SG is preferred over the OG; its ability to reproduce the light curve is seen clearly by comparing the red curves in panels (a) and (c) of Figure 7. The reverse is true under the assumption of off-peak emission from a wind nebula, shown by the blue curves. Under the assumption of magnetospheric off-pulse emission, the SG misses some emission in the wings and has too high a background level; this effect is magnified when we assume instead a PWN origin of this emission. Such details may be modified with a more physical emission model, for example by including azimuthal asymmetry in the accelerating electric field from offset polar caps, which leads to a decrease in the off-peak emission (Harding & Muslimov 2011). The OG characteristically does not reproduce the wings or higher emission in the off-peak phases when the background is derived assuming pulsar emission in the off-peak; by definition it matches the background well when the background is instead found assuming the emission is not from the pulsar itself.
Regardless of the background counts used, the best-fit values of α, and ζ for the OG and SG models are far apart, with |ζ −α| >60 degrees. The radio beam would need a width of that order in order for the radio pulse to be detectable. For the parameters of PSR J0007+7303 a model estimate of beam width is < 10 degrees (Story et al. 2007). Thus the preference for OG and SG model fits over that of TPC, along with the fitted parameter values, self-consistently offers a satisfactory explanation for non-detection of any radio pulse. A deep search using GBT (Halpern et al. 2004) yielded an upper limit that remains the lowest limiting flux for any pulsar position. The upper limit on the pseudo-luminosity of L 1400 ≃ 0.02 mJy kpc 2 is lower than the lowest pseudo-luminosity measured for a radio pulsar (L 1400 ≃ 0.035 mJy kpc 2 for PSR J1907+0602 (Abdo et al. 2010d)) and is thus restrictive enough to support the radio-quiet designation. Conversely, if one were to adopt the TPC model for this source then the nearly-equal values of α and ζ would suggest that a radio pulse should have been found.
X-ray pulsations were recently detected from the CTA 1 pulsar using the same LAT ephemeris used in the timing and light curve analysis of this paper (Lin et al. 2010;Caraveo et al. 2010). We did not include any X-ray information in our fits, but can consider whether or not our results make sense in the multiwavelength picture. The peak in the X-ray light curve shows a significant thermal component, suggestive of a hot spot on the surface. The authors estimate the size of the hot spot to be ∼ 100 m in radius, larger than the polar cap in the case of a dipole model for the pulsar; it is not clear from the data where the hot spot is located. The X-ray peak occurs at φ ∼ 0.25-0.3 (∼ 90 • -110 • ) prior to the first peak in the γ-ray light curve. Our HE models predict a similar location of the magnetic pole: the first γ-ray peak is at 90 • , and in all our model fits the pole lies 36 • -94 • before this peak; taking only the SG and best OG fits places the pole 94 • -88 • . Our best-fit models are therefore consistent with the thermal X-ray emission originating near the magnetic pole.

Thermal Component in X-rays and Neutron Star Cooling
PSR J0007+7303 is among the neutron stars (NS) where the theory of cooling confronts observation. It held that distinction before launch of Fermi, based on the candidate NS identified at that time (Halpern et al. 2004), and it continues to be treated as an exceptional case in the cooling literature (Page et al. 2009). Those references cite additional literature on heat loss models which describe how observable surface temperature or thermal flux declines with age. X-ray or extreme ultra violet surface emission in young, hot NS provides the observational test. The NS age and surface temperature must be known or constrained, that is, a limit may be useful if corrections can be characterized at least as to sign if not magnitude. Only comparatively young NS within a few kpc provide useful constraints. The only pulsar both younger and closer than PSR J0007+7303 is the Vela pulsar. Before Fermi, the age of PSR J0007+7303 was equated to the estimated CTA1 SNR age and thermal flux was estimated from the X-ray emission of the identified point-like NS candidate (Halpern et al. 2004). From this it was understood that the putative pulsar in CTA1 was cool for its age, in comparison with both theory and other young NS. This led to theoretical effort to understand why this particular NS had cooled rapidly.
Fermi results impact context information that enters into the analysis of the cooling question. The pulsar ephemeris affects the X-ray analysis used to extract the thermal flux. The Fermi results have not greatly altered best-estimates of key quantities but provide a new network of interlocking constraints that give greater robustness of the observational context information, as follows. (1) The spin period P and its derivativeṖ give spindown age of ∼13,900 yr for the pulsar. This age estimate is well within the broader range of historical estimates for the SNR age, 5,000-15,000 yr, and nearly identical to one pre-launch estimate of 1.3 × 10 4 d 1.4 yr (Slane et al. 2004). The formal propagated error in the spindown age estimate is negligible. The principal issue is validity of the dipole approximation formula but it clearly reinforces the age estimates derived other ways. (2) The same observations (P , P ) also yield spindown energy loss independently of distance. There are no inconsistencies between this energy budget and that of the pulsar plus SNR, using the accepted distance, D, of 1.4 ± 0.3 kpc, hence the pre-Fermi distance estimate used is not in conflict with new information and the distance scale factor in the age estimate just cited is consistent with unity. (3) The pulsar ephemeris has now been used to detect X-ray pulsations using XMM (Lin et al. 2010;Caraveo et al. 2010). The X-ray spectrum can be divided into pulsed and unpulsed components. While total X-ray flux is consistent with levels measured in pre-Fermi work, the part now assigned to a DC thermal component emitted from the surface is reduced by at least a factor that equates to the unpulsed fraction. Since D is unchanged, the revised bound on NS surface temperature remains as in earlier estimates or perhaps even slightly reduced. One may consider modifications due to the interstellar absorbing column and to the neutron star atmosphere causing the surface flux to depart from blackbody. These considerations properly enter into many other NS used for comparison with cooling and it is beyond scope of this paper to reanalyze the entire X-ray discussion. It is left to future research but provisionally the pre-launch thermal luminosity stands. (4) In this same connection a further point to note is that the off-pulse γ-ray flux, discovered and quantified in §2.2 of the present paper, cannot be thermal flux from the neutron star and must be taken as magnetospheric or PWN emission. In either case, that same magnetosphere or PWN component could extend downward in energy to X-rays, and in principle could provide further downward adjustment to the thermal X-ray flux from the star. This is left as an adjustment of unknown magnitude but known sign; it can only require the star to cool faster than the estimate obtained by neglecting it. (5) From P andṖ one also obtains an estimate of the stellar dipole field, B ∼ 10 13 G. PSR J0007+7303 therefore falls among the most highly magnetic neutron stars, magnetars excepted. Cooling scenarios whereby a strong magnetic field could modify the cooling curve have been described in earlier literature (for a summary see for example Yakovlev et al. (2001)). The high value of B now established for this pulsar means mechanisms whereby high B enhances cooling may merit further attention. (6) The previous section shows how model fits for OG and SG models favor small α and large ζ, and could explain the absence of radio pulsations.
The first four items in the list mean that if PSR J0007+7303 was an outlier relative to models before Fermi and relative to pre-Fermi theoretical understanding, it has become slightly more egregious relative to that prior theoretical understanding. However theory has also advanced, largely from observation of cooling in the central star in Cas A (Page et al. 2009). Data from that source are now the most definitive and constraining of any unrecycled NS. Fitting Cas A has given prominence to NS cooling models involving Cooper pairing contributions. Points (5) and (6) in the list regarding magnetic field strength and geometry may provide guidance to theory. Sufficiently strong magnetic fields can affect cooling (Yakovlev et al. 2001) and might affect heat flow in the star. Also the particular geometry in PSR J0007+7303 with low α and high ζ that operates against radio pulse detection could also lead to a misleadingly low thermal X-ray flux if the strong field conveyed internal heat flow preferentially to the magnetic polar regions. Then, a Lambertian emission pattern from the hotter poles would be anisotropically beamed away from an observer at high ζ. Further X-ray observations, combined with multi-wavelength analysis applied and field geometry modeling may shed further light on the thermal luminosity. This could be undertaken comparatively with other young, nearby pulsars such as the Dragonfly pulsar, PSR J2021+3651, which has estimated age 17 kyr, distance 2.1 kpc, dipole field 3 × 10 12 G, and has radio pulses. Such comparisons might bring out the role of the magnetic field strength and geometry in NS cooling.

Summary
PSR J0007+7303 is among the brightest γ-ray pulsars (Abdo et al. 2010c). It also is part of an interesting PWN and SNR complex that is still young. We have exploited greatlyimproved cumulative statistics from two years of Fermi -LAT data to investigate questions that can be pursued only on brighter γ-ray pulsars. Interesting aspects of the source include its having had a major glitch during the Fermi observing period, its being among the most strongly magnetic of neutron stars that are not magnetars, and its status as a prime testbed for neutron star cooling theory, a topic now receiving renewed attention because of the cooling observed in X-rays in Cas A.
Pulsar phase dependence of the light curve and spectrum have been investigated and the source has been compared both with other well-studied pulsars (notably PSR J1709-4429) and also with standard magnetospheric geometry models, revealing a clear preference for models where emission occurs high in the magnetosphere; particularly the SG model. Glitch parameters have been extracted. We have conducted a systematic search for long-term variability in the system, with negative results. Neither a change associated with the glitch nor flaring such as has recently been seen in the younger Crab Nebula (Abdo et al. 2011) has been found. However, off-pulse emission has finally been detected at high confidence and there is evidence that it is extended; hence there is now a new γ-ray component in the overall source, potentially a PWN although the possibility that it originates inside the magnetosphere is not strongly excluded. The variability of that off-pulse source is a subject to be pursued as Fermi continues to accumulate data. We have described how the parameters emerging from the γ-ray analysis (ephemeris, spindown energy loss, age, distance, and magnetic field) and from follow-on X-ray studies affect understanding of the cooling history, reinforcing the conclusion that the surface of this NS is cool for its age. This now needs to be followed up with additional X-ray analysis to further constrain the surface X-ray luminosity.
We thank Mallory Roberts and Tyrel Johnson for helpful contributions.
The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariatà l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.
Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.
This work was performed under contract with the naval research laboratory, contract N000173-08-2-C004 and was sponsored under a grant by NASA. Black cross marks the location of the pulsar. ROSAT X-ray contours are shown in black (Seward et al. 1995). Radio contours are shown in green (Pineault et al. 1997).  Table 2) a power law with exponential cutoff has been assumed. The dashed histogram represents the LAT light curve above 0.1 GeV. -The best fitting model light curves for each model, red indicating 195 counts/bin and green 246 counts/bin as the background level, and the reduced χ 2 r = χ 2 /N dof contours in α-ζ parameter space for the case of a background of 195 counts/bin. (a) Best fitting outer gap light curve (red/blue) superposed on the LAT light curve of the PSR J0007+7303 pulsar (black); the best fit parameters are given in the text. The horizontal dashed lines represent the estimated background levels, while the vertical dotted lines mark the location of the magnetic pole in the context of the outer gap model for the listed sets of parameters. (b) χ 2 r contours for outer gap light curves in which w and r have been fixed at the best-fit parameters. The filled contours show χ 2 in increments of 20 for all α and ζ. The best fit is marked by the pink triangle. The 3 σ confidence intervals for each parameter are given in Table 4. (c) Same as (a), for the slot gap model. (d ) Same as (b), for the slot gap model.  Table 2: Phase interval definitions and corresponding spectral parameters obtained from fitting the spectrum with a power law with exponential cutoff. The flux in the third column is normalized to the width of the phase bin. The systematic uncertainties are in agreement with the ones evaluated for the phase-averaged analysis. 3.60 ± 0.11 ± 0.27 4.09 ± 0.08 ± 0.31 Table 3: Phase-averaged spectral parameter and flux for PSR J0007+7303 for the two epochs around the glitch. Flux is given in 10 −  Table 4: Best-fit model parameters of the OG and SG geometrical models to the LAT light curve of PSR J0007+7303, where the background used for the fitting was 195 counts/bin for columns labeled "1" and 246 counts/bin for those labeled "2". The asymmetric error bars give the 3σ confidence intervals, derived using the scaled χ 2 as described in §3.1. There were two regions in parameter space that fell within 3σ of the best OG1 and SG2 fits; the first set of values given in the column are the best of the two regions. The quoted reduced χ 2 and f Ω correspond to the absolute best fits in each region. The parameters shown here are used in Figure 7 for the simulated light curves and χ 2 contours.