FERMI-LAT OBSERVATIONS OF THE GEMINGA PULSAR

The strong energy dependence of the PSF imposes energy-dependent ROIs that optimize the signal-to-noise ratio. Following a procedure similar to that used for the Fermi-LAT pulsar catalog paper (Abdo et al. 2009i), to study the pulse profiles we selected photons within an angle θ < max[1.6–3log₁₀(E_GeV), 1.3] deg from Geminga. Such selection provides clean light curves by limiting acceptance of the softer Galactic background.

We used the Fermi tool gtpphase to correct photon arrival times to the solar system barycenter using the JPL DE405 solar system ephemeris (Standish 1998) and to assign a rotational phase to each photon using the timing solution described in Section 3.

Figure 2 shows the light curve of Geminga above 0.1 GeV obtained with the energy-dependent cut. In order to better show the fine structure, we plot the pulse profile using variable-width phase bins, each one containing 400 events. The photon flux in each phase interval thus has a 1σ Poisson statistical error of 5%. The dashed line represents the contribution of the diffuse background, estimated by selecting photons in the phase interval ϕ = 0.9–1.0 from an annular region between 2° and 3° from the source rescaled for the solid angle and also taking into account the energy-dependent selection adopted. The light curve contains 61219 ± 284 pulsed photons and 9821 ± 99 background photons.

The pulse profile shows two clear peaks at ϕ = 0.141 ± 0.002 (P1) and ϕ = 0.638 ± 0.003 (P2). In order to reveal possible asymmetries in the peaks, we started by fitting the sharp peaks with two half-Lorentzian profiles with different widths for the trailing and the leading edge. We have chosen this function because it has a simple parameterization and appears to fit well the pulse profile of the gamma-ray light curves. We found that Geminga peaks show no asymmetries, and P1 is broader (FWHM of 0.072 ± 0.002) than P2 (FWHM 0.061 ± 0.001). We also checked if the peaks can be better fitted by a Gaussian profile, finding comparable results (P1 FWHM of 0.071 ± 0.002) and (P1 FWHM of 0.063 ± 0.001), though we cannot distinguish between a Lorentzian or a Gaussian profile. The smallest features in the pulse profile appear on a scale of 260 μs, presumably artifacts of the timing model residuals. Figure 2 also contains insets (binned to 0.00125 in phase) centered on the two peaks and on the phase interval ϕ = 0.9–1.0. This off-peak, or "second interpeak," region contains 789 ± 28 pulsed photons above the estimated background (∼1.3 × 10⁻² of the pulsed flux). This corresponds to a signal-to-noise ratio of 19σ, indicating that the pulsar emission also extends in the off-peak, as will be investigated further in Section 4.3.

Figure 3 shows the pulse profile in five energy ranges (0.1–0.3 GeV, 0.3–1 GeV, 1–3 GeV, 3–10 GeV, and >10 GeV). There is a clear evolution of the light-curve shape with energy: P1 becomes weaker with increasing energy, while P2 is still detectable at high energies. Significant pulsations from P2 are detectable at energies beyond _max ∼ 18 GeV, chosen as the maximum energy beyond which a χ² periodicity test still attains 6σ significance. We detect 16 photons above 18 GeV, not necessarily coming from the pulsar itself. No particular features appear at high energies in the bridge region between P1 and P2 ("first interpeak").

Zoom In Zoom Out Reset image size

Figure 4 shows the evolution of the P1/P2 ratio as a function of energy, plotted using variable-width energy bins. The curve depends very weakly on the bin choice; Figure 4 was made using 10,000 events per bin. A clear decreasing trend is visible, as observed in the Crab, Vela, and PSR B1951+32 γ-ray pulsars by EGRET (Thompson 2004) and now confirmed for the Vela (Abdo et al. 2009a) and the Crab pulsars (Abdo et al. 2010c) by Fermi-LAT. Adopting the same variable-width energy bins, we fit the peaks in each energy range with a Lorentz function to determine the peak center and width. Figure 5 shows the energy evolution of the FWHM of P1 and P2: both peaks narrow with increasing energy. The decreasing trend in pulse width of P1 and P2 is nearly identical. P1 has an FWHM decreasing from δϕ = 0.098 ± 0.004 to δϕ = 0.053 ± 0.008, while FWHM of P2 changes from δϕ = 0.092 ± 0.004 to δϕ = 0.044 ± 0.004 at energies greater than 3 GeV. The decrease in width with energy does not depend on the shape used to fit the peaks. Figure 8 was made using the Lorentzian fits, preferred in general because they are sensitive to asymmetric pulses. While the "first interpeak" emission is significantly detected up to 10 GeV, emission in the "second interpeak" region (between 0.9 and 1.0), not detected before, is clearly present at low energies but vanishes above ∼2 GeV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Spectral analysis was performed using the maximum-likelihood estimator gtlike included in the standard Fermi Science Tools provided by the FSSC. The fit was performed using a region of the sky with a radius of 15° around the pulsar position, selecting energies between 0.1 and 100 GeV.

We included in the fit a model accounting for the diffuse emission as well as for the nearby γ-ray sources. We modeled the diffuse foreground, including Galactic interstellar emission (Casandjian & Grenier 2008; Strong et al. 2004a, 2004b), extragalactic γ-ray emission, and residual CR background, using the models⁶¹ gll_iem_v02 for the Galactic part and isotropic_iem_v02 for the isotropic one.

In the fit procedure, we fixed the spectral parameters of all the sources between 15° and 20° from Geminga, and left free the normalization factor of all the sources within 15°. All the non-pulsar sources have been modeled with a power law as reported in the Fermi Bright Source List (Abdo et al. 2009i), while all the pulsars have been described by a power law with exponential cutoff according to the data reported in the Fermi-LAT pulsar catalog (Abdo et al. 2009i).

We integrated the phase-averaged spectrum to obtain the energy flux. The unbinned gtlike fit is described by a power law with exponential cutoff in the form:

$$\begin{equation} \frac{\it dN}{\it dE} = N_{0} E^{-\Gamma } \exp \left(- \frac{E}{E_{0}}\right) \mbox{cm$^{-2}\; $s$^{-1}\; $GeV$^{-1}$}, \end{equation} \tag{ 1 }$$

where N₀ = (1.189 ± 0.013 ± 0.070) × 10⁻⁹ cm⁻² s⁻¹ GeV⁻¹, Γ = (1.30 ± 0.01 ± 0.04), and E₀ = (2.46 ± 0.04 ± 0.17) GeV. The first uncertainties are statistical values for the fit parameters, while the second ones are systematic uncertainties. Systematics are mainly based on uncertainties on the LAT effective area derived from the on-orbit estimations, and are ⩽5% near 1 GeV, 10% below 0.1 GeV, and 20% above 10 GeV. We therefore propagate these uncertainties using modified effective areas bracketing the nominal ones (P6_v3_diffuse).

For this fit, over the range 0.1–100 GeV, we obtained an integral photon flux of (4.14 ± 0.02 ± 0.32) × 10⁻⁶ cm⁻² s⁻¹ and a corresponding energy flux of (4.11 ± 0.02 ± 0.27) × 10⁻⁹ erg cm⁻² s⁻¹.

We studied alternative spectral shapes beginning with the cutoff function exp [ − (E/E₀)^b]. The 46 gamma-ray pulsars discussed in Abdo et al. (2010a) are generally well described by a simple exponential cutoff, b = 1, a shape predicted by outer magnetosphere emission models (see Section 5). Models where gamma-ray emission occurs closer to the neutron star can have sharper "super-exponential" cutoffs, e.g., b = 2. Leaving free the exponential index b, we obtained N₀= (1.59 ± 0.13 ± 0.09) × 10⁻⁹ cm⁻² s⁻¹ GeV⁻¹, Γ = (1.18 ± 0.03 ± 0.04), E₀ = 1.58 ± 0.19 ± 0.11) GeV, and b = (0.81 ± 0.03 ± 0.06). As previously reported for the analysis of the Vela pulsar (Abdo et al. 2010b), b < 1 can be interpreted by a blend of b = 1 spectra with different cutoff energies. Figure 6 shows the results of the phase-averaged spectrum in case of b free (dashed line) and b fixed to 1 (solid line). Using the likelihood ratio test, we found that the hypothesis of b = 2 can be excluded since the likelihood of this fit being a good representation of the data is much greater than for a power-law fit (logarithm of the likelihood ratio being 396). We have also tried different spectral shapes, like a broken power law, but the fit quality does not improve (the logarithm of the likelihood ratio is 212).

Zoom In Zoom Out Reset image size

We divided the pulse profile in variable-width phase bins, each one containing 2000 photons according to the energy-dependent cut defined in Section 4.1. This choice of binning provides a reasonable compromise between the number of photons needed to perform a spectral fit and the length of the phase intervals that should be short enough to sample fine details on the light curve, while remaining comfortably larger than the rms of the timing solution (Section 3). We have performed a maximum-likelihood spectral analysis, similar to the phase-averaged one, in each phase bin assuming a power law with exponential cutoff describing the spectral shape. Using the likelihood ratio test, we checked that we can reject the power law at a significance level greater than 5σ in each phase interval. Following the results on phase-averaged analysis of Geminga, we have modeled the spectrum in each phase interval with a power law with exponential cutoff. Such a model yields a robust fit with a logarithm of the likelihood ratio greater than 430 in each phase interval. Figure 7 (below) shows the evolution of the spectral parameters across Geminga'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 3 summarizes the results of the spectral fit in each phase bin. In this case, we have fixed all 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.

Zoom In Zoom Out Reset image size

To obtain Fermi-LAT spectral points, we divided our sample into logarithmically spaced energy bins (four bins per decade starting from 100 MeV) and then applied the maximum-likelihood method in each bin. For each energy bin, we have used a model with all the nearby sources as well as Geminga described by a power law with a fixed spectral index. We have considered only energy bins in which the source significance was greater than 3σ. From the fit results, we then evaluated the integral flux in each energy bin. This method does not take energy dispersion into account and correlations among the energy bins. To obtain the points of the spectral energy distributions (SEDs), we multiplied each bin by the mean energy value of the bin taking into account the spectral function obtained by the overall fit. Figures 9–12 in the Appendix show the SEDs obtained in each phase interval. The fluxes in Y-axis are not normalized to the phase bin width, whereas in Table 3 of the Appendix the fluxes are normalized. Figure 7 shows the phase evolution of the spectral index and cutoff energy, respectively. The spectral index reaches a local minimum around P1 (ϕ ∼ 0.14–0.15) and, after a sudden increase, begins to decrease again in the "first interpeak" region, reaching a minimum of Γ ∼ 1.1 around the leading edge of P2 (ϕ ∼ 0.60–0.61). It then starts to rise again in the phase interval from P2 to the "second interpeak" region (ϕ = 0.9–1.0).

The cutoff energy evolves quite differently as a function of the rotational phase. It closely follows the pulse profile, thus confirming the observations performed by EGRET (Fierro et al. 1998), which unveiled a correlation between hardness ratio and pulse profile. As shown in EGRET data and recently confirmed by AGILE (Pellizzoni et al. 2009), the hardest component is P2: our phase-resolved scan points to a cutoff around 3 GeV and a spectral index of ∼1.0 that become softer through the peak. P1 appears to be softer, with a cutoff energy slightly greater than 2 GeV and a spectral index Γ ∼ 1.2.

The phase-resolved spectra show that Geminga's emission in the bridge (or "first interpeak") phase interval (ϕ = 0.2–0.52) is quite different from the Crab (Abdo et al. 2010c) or Vela pulsars (Fierro et al. 1998; Abdo et al. 2009a). For the Crab pulsar, the bridge emission shows no evolution and drops to an intensity level comparable to the off pulse emission, while for the Vela pulsar it varies substantially but is always seen at high energies. The "first interpeak" of Geminga, instead, becomes harder and remains quite strong at high energies, as can also be seen in Figure 3. Another difference with respect to the Vela pulsar is that Geminga does not have a third peak like the one observed at GeV energies in the Vela pulsar (Abdo et al. 2009a).

The analysis of the "second interpeak" region around ϕ = 0.9–1.0 shows significant emission up to ∼2 GeV (Figure 3). Moreover, the spectrum in this phase interval has been fit with a power law with exponential cutoff, obtaining a spectral index Γ = (1.48 ± 0.17) and E₀ = (0.87 ± 0.19) GeV, with systematic uncertainties in agreement with those evaluated in the phase-averaged analysis. A pure power-law fit can be rejected with an ∼8σ confidence level, thus confirming the presence of the cutoff. The presence of the "second interpeak" component is also visible in the maps of Figure 8, where the emission in this phase region is not visible at high energies, as expected owing to the spectral cutoff.

Zoom In Zoom Out Reset image size

Analyzing the phase evolution of the spectral parameters in Figure 7, it seems that no abrupt changes occur in this phase interval and that this emission may be related to the wings of the peaks. This fact, together with the newly detected off-peak emission, favors a pulsar origin of such "second interpeak" emission, rather than an origin in a surrounding region. The detection of off-peak emission, rendered possible by the outstanding Fermi statistics, is a novelty of Geminga's high-energy behavior.

The unprecedented photon statistics collected by Fermi-LAT allows for tighter observational constraints on emission models. The absence of radio emission characterizing Geminga clearly favors models where the high-energy emission occurs in the outer magnetosphere of the pulsar.

Polar Cap (PC) models, where the high-energy emission is located near the neutron star surface (Daugherty & Harding 1996), are unlikely to explain the Geminga pulsar, since the line of sight is necessarily close to the magnetic axis for such models where one expects to see radio emission.

The current evidence against low-altitude emission in γ-ray pulsars (Abdo et al. 2009i) can also be supplemented by constraints on a separate physical origin. In PC models, γ-rays created near the neutron star surface interact with the high magnetic fields of the pulsar, producing sharp cutoffs in the few to ∼10 GeV energy regime. Moreover, the maximum observed energy of the pulsed photons observed must lie below the γ–B pair production mechanism threshold, providing a lower bound to the altitude of the γ-ray emission. According to Baring (2004), the lower limit for the altitude of the production region r could be estimated taking advantage of the maximum energy detected for pulsed photons _max as $r\ge \left(\epsilon _{\rm max}B_{12}/1.76\; \rm GeV\right)^{\frac{2}{7}}P^{-\frac{1}{7}}R_{*}$, where P is the spin period, R_* is the stellar radius, and B₁₂ is the surface magnetic field in units of $10^{12}\; \rm G$. For pulsed photons of _max ∼ 18 GeV, we obtain r_min ⩾ 2.7R_*, a value clearly precluding emission very near the stellar surface, adding to the advocacy for a slot gap (SG) or outer gap (OG) acceleration locale for the emission in this pulsar.

OG models (Cheng et al. 1986; Romani 1996; Zhang & Cheng 2001), where the high-energy emission extends between the null charge surface and the light cylinder, the two-pole caustic (TPC) models (Dyks & Rudak 2003) associated with SG (Muslimov & Harding 2004), where the emission is located along the last open field lines between the neutron star surface and the light cylinder, or a striped wind model (Pétri 2009), where the emission originates outside the light cylinder, could produce the observed light curve and spectrum. Nevertheless, the observed peak separation of 0.5 is unlikely for a middle-aged pulsar like Geminga in the OG model, if it is true that emission moves to field lines closer to the magnetic axis as pulsars age. For the OG model, this drift leads to <0.5 peak separations. For the TPC models, 0.5 peak separation can occur in spite of this shift, that is, for all ages and spin-down luminosities.

Following the Atlas of γ-ray light curves compiled by Watters et al. (2009), we can use Geminga's light curve to estimate, for each model, the star's emission parameters, namely, the Earth viewing angle ζ_E with respect to the neutron star spin axis and the inclination angle α between the star's magnetic and rotation axes. Table 2 summarizes the observed parameters and gives the estimated beaming correction factor f_Ω(α, ζ_E), which is model-sensitive. It is given by Watters et al. (2009) as

$$\begin{equation} f_{\Omega }(\alpha,\zeta _{E}) = \frac{\int F_{\gamma }(\alpha ;\zeta,\phi)\sin (\zeta)d\zeta d\phi }{2 \int F_{\gamma }(\alpha ;\zeta _{E},\phi) d\phi }, \end{equation} \tag{ 2 }$$

where F_γ(α; ζ, ϕ) is the radiated flux as a function of the viewing angle ζ and the pulsar phase ϕ. In this equation, the numerator is the total emission over the full sky, and the denominator is the expected phase-averaged flux for the light curve seen from Earth.

The total luminosity radiated by the pulsar is then given by L_γ = 4πf_ΩF_obsD², where F_obs is the observed phase-averaged energy flux over 100 MeV and D = 250⁻⁶²₊₁₂₀ pc is the pulsar distance (Faherty et al. 2007). The estimated averaged luminosity is then L_γ = 3.1 × 10³⁴f_Ω erg s⁻¹, yielding a γ-ray efficiency $\eta _{\gamma }=\frac{L_\gamma }{\dot{E}}$ = 0.15f_Ω (d/100 pc)².

Ideally, geometrical values in Table 2 should be compared with independent estimates, coming, e.g., from radio polarization or from the geometry of the pulsar wind nebula (Ng & Romani 2004, 2008).

Owing to the lack of radio emission, the only geometrical constraints available for Geminga come from the X-ray observations which have unveiled a faint bow shock structure due to the pulsar motion in the interstellar medium (Caraveo et al. 2003) and an inner tail structure (De Luca et al. 2006; Pavlov et al. 2006), while phase-resolved spectroscopy yielded a glimpse of the geometry of the emitting regions as the neutron star rotates (Caraveo et al. 2004).

The shape of the bow shock feature constrains its inclination to be less than 30° with respect to the plane of the sky. Since such a feature is driven by the neutron star proper motion, the constraint also applies to the pulsar proper motion vector and thus, presumably, to its rotation axis, as is the case for the Vela pulsar (Caraveo et al. 2001), pointing to an Earth viewing angle ranging from 60° to 90°.

Analyzing the pulsar spectral components along its rotational phase, Caraveo et al. (2004) concluded that the observed behavior could be explained in the frame of an almost aligned rotator seen at high inclination.

However rough, such constraints would definitely favor the OG model pointing to a beaming factor of 0.1–0.15. Such a value turns out to also be in agreement with the heuristic luminosity law $\eta \simeq \left(\dot{E}/10^{33}\right)^{-0.5}$ given by Arons (1996) and Watters et al. (2009), which for the Geminga parameters should yield a value of ∼17%. For the nominal parallax distance of 250 pc, a beaming factor of 0.15 would yield a luminosity of L_γ = 4.6 × 10³³ erg s⁻¹.

We note that the TPC models, characterized by higher efficiency, would yield higher luminosity which would account for the entire rotational energy loss for a distance of ∼300 pc, well within the distance uncertainty. On the other hand, a 100% efficiency would translate into a distance of 730 pc for the OG model, providing a firm limit on the maximum source distance.

The power law with exponential cutoff describes only approximately the phase-averaged spectrum of Geminga, since several spectral components contribute at different rotational phases. The phase-resolved analysis that we have performed is thus a powerful tool for probing the emission of the Geminga pulsar.

Figure 7 shows a sudden change in the spectral index around each peak maximum. The spectrum appears to be very hard in the "first interpeak" region between P1 and P2, with an index close to Γ ∼ 1.1 and quickly softens after the peak maximum and in the "second interpeak" to Γ ∼ 1.5. Caustic models such as OG and TPC predict such behavior as a result of the change in emission altitude with energy. Sudden changes in the energy cutoff are also predicted, as is also seen for Geminga. Large variations in the spectral index and energy cutoff as a function of the pulsar phase have already been seen in other pulsars, such as the Crab pulsar (Abdo et al. 2010c) or PSR J2021+3651 (Abdo et al. 2009e).

The persistence of an energy cutoff in the "second interpeak" region suggests pulsar emission extending over the whole rotation, further supporting the TPC model for Geminga. A similar "second interpeak" has also been observed by Fermi-LAT in PSR J1836+5925, known as the "next Geminga" (Halpern et al. 2007). Although Geminga is significantly younger, the two pulsars share other interesting features, including very similar spectral indices and energy cutoffs in the phase-averaged spectrum, and comparable X-ray spectra (Abdo et al. 2010d).