FERMI LARGE AREA TELESCOPE OBSERVATIONS OF THE CRAB PULSAR AND NEBULA

The selected γ-rays were phase-folded using the timing solution described in Section 2. Photons with an angle θ < max(6.68 − 1.76log₁₀(E_MeV), 1.3)° of the radio pulsar position, R.A. = 8363322, decl. = 2201446 (J2000), are selected. This choice takes into account the instrument performance and maximizes the signal-to-noise ratio. At high energies, the background is relatively faint compared to the Crab emission, so that a radius larger than the point-spread function (PSF) can be kept.

Using this energy-dependent region, the γ-ray light curve above 100 MeV is presented in Figures 1 and 2(g). We have 22,601 γ-rays among which we estimate 14,563 ± 240 pulsed photons after background subtraction. Phase histograms in radio (from the Nançay radio telescope), optical (Oosterbroek et al. 2008), X-rays (Rots et al. 2004), hard X-rays (Mineo et al. 2006), γ-rays (CGRO COMPTEL and EGRET; Kuiper et al. 2001), and very high energy (VHE) γ-rays (MAGIC; Aliu et al. 2008) are also plotted in Figure 2. We did not search for any correlation between giant pulses and γ-ray photons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The phase 0 is taken at the maximum of the main radio peak observed at 1.4 GHz, as seen in Figures 1 and 2(a). Considering all events between 100 MeV and 300 GeV, two clear peaks P1 and P2 can be seen at phases ϕ₁ = 0.9915 ± 0.0005 and ϕ₂ = 0.3894 ± 0.0022, respectively. Hence, the peaks are separated by δϕ = 0.398 ± 0.003. P1 and P2 are asymmetric. Their shapes can be well modeled by two half-Lorentzian functions (with different widths for the leading and trailing sides). The first peak presents rising and falling edges of half-widths 0.045 ± 0.002 and 0.023 ± 0.001, respectively. P2 shows a slow rise and a steeper fall. The rising and falling edges of P2 have Lorentzian half-widths of 0.115 ± 0.015 and 0.045 ± 0.008, respectively. Hence, the γ-ray first peak leads the radio main pulse by phase 0.0085 ± 0.0005, as shown in Figure 1, where the radio profile (red line) is overlaid for comparison.

The second γ-ray peak leads the second 1.4 GHz radio pulse (interpulse) by 0.0143 ± 0.0022 in phase. The peak separation is slightly wider at 1.4 GHz than in γ-rays.

An error in these γ-radio delays can also arise from the measurement of the dispersion measure and its derivative. Following Manchester & Taylor (1977), the error on the dispersion delay in the propagation of a signal at a frequency f through the interstellar medium is

$$\begin{eqnarray} &&\Delta (\Delta t) = - \frac{\Delta {\rm DM}}{K f^2}, \end{eqnarray} \tag{ 1 }$$

where ΔDM takes into account the error on the measurement of DM and its derivative, and K = 2.410 × 10⁻⁴ MHz⁻² cm⁻³ pc s⁻¹ is the dispersion constant. This yields a formal uncertainty of 1.4 μs, which is significantly smaller than the 21.1 μs accuracy of the overall timing solution, and therefore leads to an error of 0.0006 in phase on the γ-radio delay.

The presence of a radio feature referred to as low frequency component (LFC) by Moffett & Hankins (1996) can be noticed at phase 0.896 ± 0.001 on the radio light curve obtained at 1.4 GHz as seen in Figures 1 and 2(a). This peak is assumed to be near the closest approach of the magnetic axis. The first γ-ray peak lags the LFC by 0.095 ± 0.002 in phase.

Figure 3 shows the light curves in five energy bands, covering the 100 MeV–300 GeV interval while Table 1 reports the evolution of the positions of the peak maxima (ϕ₁ and ϕ₂ for P1 and P2, respectively) and their half-widths (HW), for the energy bins between 100 MeV and 10 GeV. The photon number counts above 10 GeV were not sufficient to fit the peak profiles. The phases of the first (P1) and second (P2) peaks do not show any significant shift with energy. Both become narrower when the energy increases, showing in particular a steepening in the P2 falling edge.

Zoom In Zoom Out Reset image size

Table 1 also presents the energy dependence of the relative weight of the two peaks. The diffuse and nebular background photon density has been first estimated in the 0.52–0.87 phase interval, then renormalized and subtracted so as to determine the number of pulsed photons in both peaks. P1 and P2 are here defined in the 0.87–1.07 and 0.27–0.47 phase intervals, respectively. As for the Vela pulsar (Abdo et al. 2009a), the ratio P1/P2 decreases with increasing energy, especially above a few GeV.

We define the off-pulse window as the 0.52–0.87 phase range, due to the bright emission of the pulsar in the rest of the phase. In the light curve above 10 GeV, we can notice an enhancement indicating a potential third peak at phase ∼ 0.74, coincident with the radio peak observed between 4.7 and 8.4 GHz and referred to as High Frequency Component 2 (HFC2) by Moffett & Hankins (1996). The excess above the background level (estimated at 2.10 counts per bin, with a bin width of 0.02 in phase) is 13.8 photons in the off-pulse interval. The statistical significance of this third peak, 2.3σ, is therefore too low to claim a definite detection and a third peak will not be considered separately in the analysis of the Crab Nebula.

Figure 4 shows the counts maps of pulsed and nebular emission in a 15° × 15° region centered on the pulsar radio position, for different energy bands. The nebular emission seen in the off-pulse window has been renormalized to the total phase (bottom row) and subtracted from the whole phase emission, to obtain the maps presenting the pulsed emission only (top row). The positions of the pulsar and nebula are coincident to within our angular resolution and the nebula appears as a point-like source. While the pulsed emission dominates in the on-pulse window, the nebula stands out in the off-pulse interval from the emission of the diffuse background at high energy only.

Zoom In Zoom Out Reset image size

The spectral analysis of the γ-ray emission of the Crab Nebula was performed using a maximum-likelihood method (Mattox et al. 1996) implemented in the Fermi Science Support Center science tools as the "gtlike" code. This fits a source model to the data along with models for the instrumental, extragalactic, and Galactic backgrounds. We used an updated instrument response function (IRF), Pass6_v3, that corrects a pileup effect identified in orbit. We selected photons in the 0.52–0.87 pulse phase window in a 20° region around the pulsar radio position. Owing to uncertainties in the instrument performance still under investigation at low energies, only events in the 100 MeV–300 GeV energy band are analyzed. The Galactic diffuse emission is modeled using GALPROP (Strong et al. 2004a, 2004b) updated to include recent H i and CO surveys, more accurate decomposition into Galactocentric rings, and many other improvements. The GALPROP run designation for our model is 54_59Xvarh7S. The instrumental background and the extragalactic radiation are described by a single isotropic component with a power-law shape. Sources nearby the Crab with a statistical significance larger than 5σ are extracted using the analysis procedure described in Abdo et al. (2009b) but with 6 months of survey data, and taken into account in the study.

The systematic errors on the spectral parameters are dominated by the uncertainties in the LAT IRFs. We bracket the energy-dependent effective area with envelopes above and below the nominal curves by linearly connecting differences of (10%, 5%, 20%) at log(E) of (2, 2.75, 4), respectively. This yields the systematic errors cited below.

In parallel to the standard analysis, we have also evaluated the spectrum using an unfolding method based on Bayes' theorem (D'Agostini 1995; Mazziotta 2009), which allows the reconstruction of the true energy spectrum from the observed one taking into account the dispersions introduced by the IRF and without assuming any model for the spectral shape. The results from this analysis are consistent with those from the likelihood analysis.

Using EGRET observations, de Jager et al. (1996) reported that the IC component dominates above ∼ 200 MeV whereas the synchrotron component is more significant at lower energies. Hence, the selected γ-ray photons should allow the study of both the fall of the synchrotron and the rise of the IC radiation.

The spectrum of the Crab Nebula between 100 MeV and 300 GeV is well described by the sum of two power-law spectra. As seen on the spectral energy distribution in Figure 5, one of the components decreases while the second increases with energy. We identify them as the falling edge of the synchrotron component and the rising edge of the IC component, respectively. The nebular spectrum can be modeled with the following function:

$$\begin{eqnarray} &&\frac{dN}{dE} = N_{\rm sync} (E_{\rm GeV})^{-\Gamma _{\rm sync}} + N_{\rm IC} (E_{\rm GeV})^{-\Gamma _{\rm IC}} \, \rm {cm^{-2} \, s^{-1} MeV^{-1}},\!\vspace*{6pt} \end{eqnarray} \tag{ 2 }$$

where N_sync = (9.1 ± 2.1 ± 0.7) × 10⁻¹³ cm⁻² s⁻¹ MeV⁻¹, N_IC = (6.4 ± 0.7 ± 0.1) × 10⁻¹² cm⁻² s⁻¹ MeV⁻¹ are the prefactors determined on 35% of the total phase, and Γ_sync = (3.99 ± 0.12 ± 0.08) and Γ_IC = (1.64 ± 0.05 ± 0.07) the spectral indices of the synchrotron and IC components. While the power-law index Γ_IC for the IC component provides a measure of the index of the mean electron/positron energy spectrum in the nebula, the synchrotron index Γ_sync possesses much less physical information, being just an indication of the steepness of the quasi-exponential turnover of the synchrotron component that peaks below the LAT energy window. Adopting a power-law fit to the synchrotron contribution apparent in the 100–400 MeV range is therefore a useful mathematical convenience. The corresponding flux above 100 MeV and renormalized to the total phase is (9.8 ± 0.7 ± 1.0) × 10⁻⁷ cm⁻² s⁻¹. The first error is statistical, whereas the second is systematic.

Figure 5 shows the spectral energy distribution in E$^2 \frac{dN}{dE}$ of the Crab Nebula renormalized to the total phase. The Fermi-LAT spectral points were obtained by dividing the 100 MeV–300 GeV range into logarithmically spaced energy bins and performing a maximum likelihood spectral analysis in each interval, assuming a power-law shape for the source. Above 5.5 GeV, the width of the energy intervals is multiplied by 3 to reduce the statistical uncertainties. These points, providing a model-independent maximum likelihood spectrum, are overlaid with the fitted model described above over the total energy range (black curve). The fit of the synchrotron (purple dashed line) and IC (blue dash-dotted line) are also represented. The spectral points and the model agree well. Statistical errors (black error bars) and the overall errors (red error bars) are plotted for the Fermi points. The EGRET spectral points are represented on the same plot. As in the case of the spectrum of the Vela pulsar (Abdo et al. 2009a), derived using an earlier set of response functions, Pass6_v1, markedly different from Pass6_v3 at low energies, the LAT spectral points at high energy indicate a lower flux in comparison to EGRET. However, it can be noticed that the Fermi flux is higher than the EGRET flux, in the low energy band dominated by synchrotron radiation.

Zoom In Zoom Out Reset image size

de Jager et al. (1996) found evidence in the EGRET data that the Crab synchrotron cut-off energy varied on timescales of the order of a year. We do not see significant variation in either the synchrotron or IC components in our more limited data span on timescales of 1, 2, or 4 months. As shown in Figure 5, a difference in flux is observed between EGRET and Fermi-LAT in the energy band dominated by synchrotron radiation as well as at higher energies (above 1 GeV). Even if variability in the synchrotron tail could be expected between EGRET and LAT, the lifetimes of the electrons producing gamma-rays via IC scattering are comparable to the remnant age, implying that the IC component should be steady in time. For these reasons, the flux change seen in the synchrotron component between EGRET and Fermi-LAT cannot be considered as significant.

The photon counts at high energy are too few for a significant cut-off or break to be seen in the flux distribution of the IC component. No cut-off or break energy can be determined at low energy for the synchrotron component using the LAT data only.

Photons from both on- and off-pulse intervals are now considered to analyze the pulsed emission. The spectral parameters of the Crab Nebula mentioned in the previous section have been renormalized to match the total phase interval and fixed to perform the spectral analysis of the Crab Pulsar.

After testing different functional forms to describe the spectrum of the pulsar, we found the best fit to be given by an exponential cut-off power-law shape:

$$\begin{eqnarray} &&\frac{dN}{dE} = N_o (E_{\rm GeV})^{-\Gamma } e^{-E/E_{c}} \, \rm {cm^{-2} \, s^{-1} MeV^{-1}}, \end{eqnarray} \tag{ 3 }$$

where N_o = (2.36 ± 0.06 ± 0.15) × 10⁻¹⁰ cm⁻² s⁻¹ MeV⁻¹ is the prefactor, Γ = (1.97 ± 0.02 ± 0.06) the spectral index, and E_c = (5.8 ± 0.5 ± 1.2) GeV the cut-off energy of the distribution. The integral flux above 100 MeV is equal to (2.09 ± 0.03 ±  0.18) × 10⁻⁶ cm⁻² s⁻¹. These results are consistent with the pulsed spectrum derived from the unfolding analysis.

Figure 6 shows the spectral energy distribution of the Crab Pulsar over the whole pulse period compared to EGRET spectral points. Results of both experiments agree well in the 100 MeV–8 GeV energy range, i.e., even at low energies, where such consistency is not observed in the spectrum of the Crab Nebula. The larger energy band covered by the LAT and its better sensitivity allows us to determine the cut-off energy of the spectrum, which was not possible with EGRET.

Zoom In Zoom Out Reset image size

We also attempted to fit the data using a power law with a generalized cut-off of the form $e^{-(E/E_c)^b}$ and found b = (0.89 ± 0.12 ± 0.28) with a likelihood value that is not significantly better than that obtained in the case of a simple exponential b = 1 cut-off. We compute the probability of incorrect rejection of other spectral shapes using the likelihood ratio test. For instance, if only statistical errors are included, the power law and hyper-exponential b = 2 hypothesis shapes are rejected at a level of 10.7σ and 4.9σ, respectively.

The large number of photons detected from the Crab allows a detailed phase-resolved spectroscopic study of its emission. Therefore, the pulse profile is divided in several intervals. The phase bins are chosen so as to contain ∼1000 pulsed photons in the energy-dependent region defined in Section 4.1. A maximum-likelihood spectral analysis is performed in each pulse phase interval, assuming a power-law and an exponential cut-off power-law shape to describe the pulsed emission.

The definition of the phase intervals is given in Table 2 along with the spectral results. The last column lists the significance of the improvement obtained when using an exponential cut-off power law instead of a pure power law, in terms of χ², if only statistical errors are included.

The corresponding spectral energy distributions are presented in Figure 7, where the horizontal error bars delimit the energy intervals. 90% C.L. upper limits were computed when the statistical significance of the energy interval was lower than 3σ. Figure 8 summarizes the phase dependence of the variation of the spectral parameters, spectral index, cut-off energy, and integral flux above 100 MeV. The vertical error bars take into account both the statistical and systematic uncertainties in the spectral parameters, while the horizontal error bars delimit the phase intervals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One observes a slight steepening of the spectrum with phase in the interval corresponding to P1 (0.955–0.098), with averaged values of spectral index and cut-off energy close to ∼ 1.9 and ∼ 3 GeV, respectively.

The 0.098–0.286 pulse phase interval presents the hardest spectrum with a spectral index of 1.49 ± 0.09 ± 0.05. This result is consistent with the spectrum of the "bridge," as defined in Fierro et al. (1998).

The spectral indices of the two peaks are the same, within the error bars, but the second peak (0.286–0.410) is characterized by a cut-off energy apparently larger than that of P1. This difference is consistent with the decrease of the P1/P2 ratio with the energy, especially above a few GeV, observed in Figure 3.

Finally, the 0.410–0.520 phase bin has the softest spectrum of the total pulse phase interval with a spectral index of 2.28 ± 0.08 ± 0.10. This explains the trend seen in the phase histograms: in Figure 3, the right edge of the second peak falls with increasing energy.

The Crab Nebula is detected across the whole electromagnetic spectrum from radio to very high energy γ-rays. The total spectral energy distribution of this source is shown in Figure 9, from soft to very-high energy γ-rays. The spectral points obtained with the LAT data analysis are also represented (red points).

Zoom In Zoom Out Reset image size

With a spectral index of Γ_IC = (1.64 ± 0.05 ± 0.07), the LAT results on the rising edge of the IC component are consistent with EGRET (1.85^+0.65_−0.5; de Jager et al. 1996). As can be observed in Figure 9, the highest part of the LAT spectrum links up satisfactorily to the lower energy bound of the Cherenkov data points. Using the LAT spectral parameters scaled to the full pulse phase, we obtain a flux at 77 GeV of (1.18 ± 0.22 ± 0.37) × 10⁻¹⁴ cm⁻² s⁻¹ MeV⁻¹, which agrees with the MAGIC differential flux at this energy of (1.14 ± 0.27 ± 0.34) × 10⁻¹⁴ cm⁻² s⁻¹ MeV⁻¹. The Cherenkov and Fermi-LAT data now cover the entire IC peak, as can be seen in Figure 9, and a break is expected at ∼ 100 GeV. Although no significant cut-off is observed in the LAT data with the current statistics, the determination of its position with an increased Fermi-LAT data sample would help the calibration of the Cherenkov telescopes, as discussed in Bastieri et al. (2005).

The IC scattering of relativistic electrons on the synchrotron, far infrared, and CMB radiation fields is considered to be the most probable mechanism for the production of γ-rays above 1 GeV. However, using a sophisticated approach carried out in the framework of the MHD flow of Kennel & Coroniti (1984), Atoyan & Aharonian (1996) have commented on the apparent deficit of GeV photons in their calculations. Taking into account both EGRET and Cherenkov results and assuming a mean magnetic field that reproduces the very high energy spectrum, they proposed that the high γ-ray flux observed by EGRET in comparison to their model is due to the enhancement of the bremsstrahlung emission from electrons captured in dense filaments. Figure 9 presents the broadband energy spectrum of the Crab Nebula together with the IC model predictions from Figure 14 of Atoyan & Aharonian (1996) for three different values of the mean magnetic field for the nebula. In view of the results obtained with the LAT, modeling the data does not require any additional emission component. The Fermi-LAT, in combination with the Cherenkov observations above 100 GeV, are in good agreement with the γ-ray flux predicted from simple IC scattering when the magnetic field lies between 100 μG and 200 μG, i.e., below the canonical equipartition field of the Crab Nebula of 300 μG. This result is consistent with the estimate of the magnetic field strength B ∼ 140 μG obtained by Horns & Aharonian (2004).

Concerning the low energy part of the nebular spectrum, the LAT spectral points, combined with COMPTEL's (taking into account statistical errors only for the latter), can be fitted with a power law with an exponential cut-off, following de Jager et al. (1996). The cut-off energy is estimated at E_c,sync = (97 ± 12) MeV. The higher value of this energy compared to that of de Jager et al. (1996) is due to the larger flux obtained with Fermi than by EGRET for the synchrotron component. The fit is represented with a blue dashed curve in Figure 9.

The high-quality statistics obtained with the Fermi-LAT both on the light curve and the spectrum of the Crab Pulsar, allow a more detailed comparison with theoretical models than previously possible. Currently, there are two classes of models that differ in the location of the emission region. The first comprises PC models which place the emission near the magnetic poles of the neutron star (Daugherty & Harding 1996). The second class consists of the OG models (Romani 1996), in which the emission extends between the null charge surface and the light cylinder, and the two-pole caustic (TPC) models (Dyks & Rudak 2003) which might be realized in SG acceleration models (Muslimov & Harding 2004), in which the emission takes place between the neutron star surface and the light cylinder along the last open field lines.

Observations of the time delay between emission at different wavelengths have been reported previously: Table 3 summarizes the delay of the radio main pulse with respect to the first peak seen from optical to high energy γ-rays. The LAT has a timing accuracy better than 1 μs (Abdo et al. 2009c) and thus enables an accurate estimation of the absolute positions of the γ-ray peaks: it was shown in Section 4.1 that the first γ-ray peak leads the radio main pulse by 0.0085 ± 0.0005 ± 0.0006 in phase, or (281 ± 12 ± 21) μs in time. Taking into account the presence of the LFC at phase 0.896 ± 0.002, the first radio peak leads the γ-ray peak by phase 0.095 ± 0.002 in phase. Observations of the evolution of the peak positions with the energy allow detailed studies of the emission regions in the magnetosphere. The delay of the radio peaks compared to other wavelengths (optical, X- and γ-rays) gives another constraint in the modeling of the emission processes taking place in pulsar magnetospheres. In particular, Kuiper et al. (2003) reported that, in the framework of the three-dimensional OG model developed by Cheng et al. (2000), one can reproduce the delay of the radio main pulse by shifting the production site of the radio emission inward toward the neutron star relative to that of high energy photons.

In the PC models, γ-rays created near the neutron star surface interact with the intense magnetic fields resulting in a sharp turnover in the few to 10 GeV energy range, while OG and SG models predict a simple exponential cut-off. Furthermore, the maximum energy of observed pulsed photons must lie below any γ–B pair production turnover threshold, providing a lower bound to the altitude of emission. We can use the observed phase-averaged cutoff energy (∼ 6 GeV) to estimate a minimum emission height as r ⩾ (_maxB₁₂/1.76 GeV)^2/7P^−1/7 R_*, where _max is the unabsorbed photon energy, P is the spin period, and the surface field is 10¹²B₁₂ G (Baring 2004). Using the parameters of the Crab Pulsar (P = 33 ms, B₁₂ = 3.78), one obtains r > 3.4 R_* which precludes emission near the stellar surface. Since we see pulsed photons up to almost _max ≈ 20 GeV in the two main pulse peaks, emission at these phases must arise at r > 4.8 R_*, with a strict lower bound of r > 3.7 R_* applying to the choice of 8 GeV, the lower energy in the highest data point window for the phase-resolved spectra in Figure 7. A similar lower bound to the emission altitude was recently reported by the MAGIC Collaboration using the hyper-exponential cutoff energy observed on the Crab Pulsar spectrum (Aliu et al. 2008). We should note here that the cut-off energy derived by the MAGIC Collaboration for a simple exponential cut-off (17.7 ± 2.8 ± 5.0) GeV is higher than the one obtained with the Fermi-LAT data, E_c = (5.8 ± 0.5 ± 1.5) GeV. However, the cut-off energy obtained with the LAT using the softer EGRET spectrum (γ = 2.022) as done by MAGIC is within the uncertainties of the MAGIC value.

To estimate the pulsed high energy γ-ray efficiency η of a pulsar, one needs to know the total luminosity radiated L_γ. It can be estimated using L_γ = 4πf_ΩF_obsD², where F_obs = (1.31 ± 0.02 ± 0.02)× 10⁻⁹ erg cm⁻² s⁻¹ is the observed phase-averaged energy flux over 100 MeV, D = (2.0 ± 0.2) kpc is the distance to the pulsar, and f_Ω is a correction factor that takes into account the beaming geometry, depending upon the magnetic inclination angle α and the Earth viewing angle ζ from the rotation axis. For the Crab Pulsar, ζ is estimated to be (63° ± 2°) from X-ray observations of the Crab Nebula torus (Ng & Romani 2008). The estimated value of α depends on the emission model used to interpret the data (Watters et al. 2009). The SG or TPC model best reproduces the observed pulse profiles for α∼ 55°–60°, whereas for the OG model α ∼ 70° gives the best result. Optical polarization measurements yield α estimates consistent with these (Slowikowska et al. 2009). The corresponding correction factors f_Ω are then equal to 1.1 for the TPC and 1.0 for the OG models. This yields a luminosity of (6.25 ± 0.15 ± 0.15) × 10³⁵ erg s⁻¹ above 100 MeV. This value is consistent with the heuristic luminosity law mentioned in Arons (1996) and Watters et al. (2009), according to which $\eta \propto \dot{E}^{-1/2}$ and verified by several γ-ray pulsars such as Vela, PSR J2021+3651 (assuming a distance of the order of 2–4 kpc), Geminga, CTA1, etc. For a neutron star moment of inertia of 10⁴⁵ g cm², the pulsed high-energy γ-ray efficiency η can be derived from the luminosity and the spin-down power $\dot{E}$: $\eta = L_{\gamma } / \dot{E} = (1.36 \pm 0.03 \pm 0.03) \times 10^{-3}$ above 100 MeV.

Knowing the value of the Earth viewing angle ζ of (63° ± 2°), the main peak separation, which is of the order of 40% of the phase, would be expected to be smaller in the PC model (Watters et al. 2009). This mismatch persists even if moderate altitudes are considered: the open field line cone opening angle enlarges to around 123 for the Crab at altitudes of around 6 stellar radii, the minimum bound inferred above from the observed absence of magnetic pair attenuation below 20 GeV. This opening angle is still somewhat too small according to the Watters et al. (2009) analysis to generate the observed main peak phase separation.