THE GALACTIC POPULATION OF YOUNG γ-RAY PULSARS

The non-recycled γ-ray pulsars almost exclusively have ages <10⁶ yr, substantially less than those of their parents, dominantly ∼10 M_☉ B stars. Thus, it is not unexpected that their distribution along the Galactic plane correlates with young star-forming regions (Yadigaroglu & Romani 1997). To exploit this correlation we need a detailed model of Galactic birth sites, namely, high-mass OB stars. On the largest scale we follow Yadigaroglu & Romani (1997) in using the free electron distribution (here we use the updated model of Cordes & Lazio 2002) as a proxy for the gas and ionizing photons associated with massive stars. This model includes the Galactic spiral arms, the 3.5 kpc ring and an exponential thin disk component with a scale height of 75 pc, which we take here as the OB star plane distribution. High-mass stars also show a "runaway" component, from disrupted binaries, giving roughly 10% of the stars peculiar velocities greater than 30 km s⁻¹ (Dray 2006). We account for this component with an exponential thick disk, as in the Cordes & Lazio model, but with a scale height of 500 pc.

The LAT pulsar population is apparently quite local (mostly ≲2 kpc) and on these scales we can employ more detailed information on the parent star distribution. In particular, the Hipparcos catalog gives a nearly complete sample of O and B stars to ∼500 pc with useful parallax information. On somewhat larger scales, catalogs of OB associations (Mel'nik & Efremov 1995) give the size and locations of massive star concentrations. We have used these data to generate a three-dimensional model of the Galactic O and B0-B2 star density. Locally, we use the direct Hipparcos sample of such stars with a ⩾2σ measurement of parallax (smoothed by a 50 pc spherical Gaussian). At 375 pc we make a smooth transition to a sample with uniform disk distributions (thin + thick runaway) augmented by the cluster OB star contribution (smoothed over 100 pc). In turn, this population merges smoothly to the large scale (spiral arm plus disk) distribution with a transition at 1.7 kpc. The overlap between the individual stars and clusters and between the clusters and the spiral arms allows us to match all components to the normalization determined by the local, complete OB star sample. The massive stars (i.e., pulsar birth sites) drawn from this distribution are shown in Figure 1 projected to the nearby Galactic plane. Although the uniform component is substantial, large OB concentrations and clusters are apparent.

Zoom In Zoom Out Reset image size

Figure 1 also shows the projected location of the nearby LAT pulsar sample. For many of these objects the distance estimates are quite poor (range shown by radial lines). Thus while several pulsars are superimposed on massive star associations, the poor distance constraints prevent confident assignment. We can, however, form a "figure of merit" to test agreement between the distribution of observed pulsar positions and positions from our detailed model simulation. This is an overlap sum over the set {i} of model birth sites with a weight determined from a Gaussian distribution about the uncertain location of each observed pulsar, j. This distribution combines the (generally large) uncertainty in the pulsar radial distances with a transverse Gaussian spread of 50 pc to account for the smoothed birth-site distribution and the offset due to pulsar proper motion. The result is

$$\begin{equation} {\rm FoM} = \sum _{j}\,\sum _{i}e^{-\frac{(l_j-l_i)^2+(b_j-b_i)^2}{2\tan ^{-1}(50\,{\rm pc}/d_j)^2}} \times e^{\frac{-(d_j-d_i)^2}{2\sigma _{d_j}^2}}. \end{equation} \tag{ 1 }$$

Errors were estimated by bootstrap analysis.

Comparisons were made between this smoothed simulation and the LAT pulsars. The detailed model gives a value of 156 ± 33. This is slightly, but not significantly, better than the values found from pulsars distributed according to the Cordes & Lazio (2002) spiral arm model (136 ± 22) or a simple uniform exponential disk alone (147 ± 26). We will thus use our detailed model in further analysis, although we caution that it is not yet demanded by the data. Improved pulsar distance estimates and/or analysis that uses observed pulsar proper motions to correct back to birth sites can substantially improve the utility of the detailed Galactic OB distribution.

One relatively robust result from the OB star modeling is a connection with Galactic core collapse rates. Summing up the local, nearly complete, Hipparcos sample for each luminosity class and using each class's nuclear evolution lifetime (Reed 2005), we can calculate the supernova rate contribution from each. We find, including stars from B2-O5, a local supernova rate of 40 SN kpc^-2 Myr⁻¹ (very high mass stars which might produce black holes contribute ∼1% of this rate). In turn our model allows us to extrapolate this to the Galaxy as a whole, yielding a Galactic O-B2 star birthrate (i.e., SN rate) of 2.4/100 yr⁻¹. This compares very well with other recent estimates, e.g., 2.30 ± 0.48/100 yr⁻¹ from nearby supernovae (Li et al. 2010), 1.9 ± 1.1/100 yr⁻¹ from ²⁶Al emission (Diehl 2006), and ∼1–2/100 yr⁻¹ from a similar analysis of nearby massive stars (Reed 2005). Thus, our model is well normalized and can be used to connect O-B2 progenitors with their γ-ray pulsar progeny.

The new-born pulsar population depends critically on the distribution of spin periods at birth. This is especially important for understanding the γ-ray sample; the bulk of the radio population is sufficiently spun down to retain little memory of this adolescent phase. We describe the birth spin periods with a single parameter truncated normal distribution:

$$\begin{equation} {\rm Prob} (P)\propto e^{-\frac{(P-P_0)^2}{P_0^2}}. \end{equation} \tag{ 2 }$$

Here, the spread equals the mean and we truncate at a minimum P = 10 ms.

Faucher-Giguere & Kaspi (2006) find a characteristic birth spin period of P₀ = 300 ms for the radio population. However, we find that this produces very few short period, γ-detectable pulsars. Indeed, one misses producing the observed γ-ray pulsar numbers by a very large factor unless the birthrate is nearly four times that allowed by the observed OB stars and supernovae (i.e., more than 10σ higher than the value estimated by Li et al. 2010). Accordingly, the birth spin period distribution, for the γ-ray pulsars at least, must extend to much lower values.

It is, however, possible to reconcile this with the radio results. We find that the requirement for a large P₀ in the radio sample is a consequence of the assumed radio luminosity law. We consider here three radio pseudo-luminosity (L = S d²) laws. The first is a power-law distribution that is independent of any other pulsar parameters (Lorimer et al. 2006). Next, Faucher-Giguere & Kaspi (2006) recommend a log-normal distribution whose central value scales as a power law with spin-down luminosity. Finally, we consider the broken power-law distribution of Stollman (1987) for which the power-law scaling of the central value stops and is flat above a cutoff (at $\dot{E} \approx 10^{34}\,{\rm erg\ s}^{-1}$). We will refer to these functions as L_R Flat, L_R Power Law, and L_R Broken Power Law, respectively.

To illustrate the effect of the initial spin and radio luminosity laws, we compare with the young pulsar populations: the LAT pulsar sample (all sky) and the high energy ($\dot{E}\ >\ 10^{33}\,{\rm erg\ s^{-1}}$) PMBS detections (within the PMBS footprint). To gauge pulsar detectability we adopt a standard radio beamwidth:

$$\begin{equation} \rho =5.8\,\deg P^{-1/2} \end{equation} \tag{ 3 }$$

(Rankin 1993), with the assumption of a circular beam of uniform (integrated) radio flux directed along the magnetic dipole axis; we discuss briefly later the possibility that young objects have even wider radio beams from high altitude (Karastergiou & Johnston 2007). We simulate the Galactic pulsar population as summarized above, assign a pseudo luminosity, check pulsar detectability at the modeled beam orientation and distance, compute the cumulative radio luminosity functions of the model detections (Figure 2), and compare with the observed samples. In all cases the "Flat" luminosity law does not reproduce the observed luminosities, so we discard this form. The broken power-law and power-law distributions are quite similar for the radio sample (upper panel), but differ for the more energetic radio-detected pulsars in the γ-ray sample. We remind the reader that we have not followed the detailed selection effects of the radio surveys. However, these are unlikely to affect these conclusions. As an example, we applied a roll-off of the sensitivity with period, as in Dewey et al. (1985). Here, S_min = S₀[9w_e/(P − w_e)]^1/2 with the effective pulse width w_e = [(0.1P)² + 10 ms]^1/2 to approximate the PMBS losses at DM ≈ 200 pc cm^-3. This caused only a small (∼6%) decrease in the number of detections for the most energetic pulsars (${\rm Log}({\dot{E}}) > 36.5$) and negligible loss (<2%) from the sample as a whole. In particular, we recover the Faucher-Giguere & Kaspi (2006) result that, with the Power Law model, P₀ = 300 ms provides a very good match to the PMBS radio luminosity function.

Zoom In Zoom Out Reset image size

However, we see that the Fermi sample (Figure 2, lower panel) shows some discrimination between the models, especially for the high-luminosity pulsars which dominate this set. The numbers of high ${\dot{E}}$ pulsars are particularly sensitive to P₀. Accordingly, we have explored this sensitivity by comparing the modeled $N(\dot{E})$ of pulsars in the energetic ($\dot{E} > 10^{33}\,\rm erg\, s^{-1}$) PMBS and Fermi samples as we vary P₀. Our statistic is a χ² comparison of the pulsar detections in half-decade bins of $\dot{E}$. We also use the additional measurement of the Galactic core-collapse rate (and error estimate) from Li et al. (2010), fitting to the total birthrate. The results are shown in Figure 3.

Zoom In Zoom Out Reset image size

As expected, this shows that with the Power Law luminosity model the radio sample prefers long initial periods, but P₀ = 300 ms is completely unacceptable for the γ-ray sample. This is because of the very small birthrate of energetic pulsars for this large P₀. One can, of course, imagine a γ-ray emissivisity law allowing high efficiency for detection of these very rare pulsars. The cost is that such pulsars are then extremely luminous and seen at very large distance. We find that the γ-ray pulsar sample gives such a sensitive probe of the high ${\dot{E}}$ population that changes sufficient to accommodate P₀ as large as 300 ms predict a pulsar sample in strong disagreement with the observed population (see further discussion in Section 4 for a specific example).

However, if we adopt the Broken Power Law radio luminosity law, we find the preferred value of P₀ for the PMBS sample decreases and that the χ² at the P₀ ≈ 50 ms demanded by the γ-ray pulsars shows only a small increase from the minimum value. We suspect that a complete treatment of the radio survey selection effects which should decrease the detectability of short P radio pulsars along the Galactic plane would further improve the agreement at small P₀.

In sum the radio and γ-ray pulsars can come from the same population if P₀ ≈ 50 ms and the extreme increase in radio luminosity for the shortest period pulsars is mitigated with the Broken Power Law. A second quantitative comparison of this agreement is shown in Figure 2, where we quote the model-data Kolmogorov–Smirnov comparison for the two pulsar samples for each luminosity law (here P₀ = 50 ms). The Broken Power Law indeed provides the best match to the detected radio luminosity function, especially for the PMBS sample. For the smaller Fermi sample the improvement is not large. We have confirmed that both this test and the $N(\dot{E})$ test are insensitive to the chosen γ-ray model (see Section 3). We thus adopt the "Broken" luminosity law as the best fit to the data.

As a consistency check, applying our adopted beaming and flux law to the population, we find that the observed number of PMBS pulsars corresponds to a full Galactic pulsar birthrate of 1.43/100 yr^-1. This is in reasonable agreement with earlier radio population estimates and with the core-collapse estimates in Section 2.1 so we may turn our attention to the γ-ray-detection criteria.

To complete the population synthesis we require a γ-ray pulsar flux and beaming model. Such models also depend on the pulsar's energetics and geometry, and a major goal of the data comparison is to constrain such dependence. All of the γ-ray models employed in this work utilize the same heuristic efficiency law

$$\begin{equation} \eta =(10^{33} \,{\rm erg\, s^{-1}}/{\dot{E}})^{1/2}, \end{equation} \tag{ 4 }$$

with γ-ray production ceasing for objects with $\dot{E} < 10^{33}\,\rm erg\,s^{-1}$. This law has theoretical motivation (Arons 2006) and provides a reasonable match to the observed pulsar efficiencies (Abdo et al. 2010a). The portion of the magnetosphere producing γ-rays and hence the fraction of the sky covered by the beam, however, differs appreciably between models.

In this paper, we treat several versions of the two popular vacuum magnetosphere beaming models. The first is the standard Outer Gap (OG) model, with emission in each hemisphere produced by field lines above a single pole. In its basic form the γ-ray production is dominated by the upper boundary of a vacuum gap which is bounded by the last closed field lines and runs from the null charge surface (Ω · B = 0) to the light cylinder R_LC = cP/2π. This gap has a spin-down-dependent thickness approximated by

$$\begin{equation} w=\eta =({10^{33}{\rm erg\ s}^{-1}}/{\dot{E}})^{1/2}, \end{equation} \tag{ 5 }$$

where w is the fractional distance from the last closed field lines to the magnetic axis, and the emission peaks ∼w away from the closed zone. As the pulsar spins down, w grows; this causes the radiating particles to approach the central field lines, the null charge crossing to move to high altitudes, and the γ-ray emission to become more tightly confined to the spin equator. To compute light curves from this simple version, we compute photon emission from particles traveling in a sheet centered on w, with a Gaussian cross section. As ${\dot{E}}$ decreases, we might expect pair production to become more difficult and this sheet to thicken. An extreme assumption is that emission arises from particles spanning the full thickness w inward from the closed zone. Models computed for this assumption are denoted full gaps (OG_FG). In both versions, this is essentially a "hollow" cone model which tends to produce two caustic peaks, of varying separation, with a bridge in between.

The other popular beaming model has emission extending from the star surface to high altitude, so that both poles are visible in both hemispheres. By truncating the emission somewhat before the light cylinder it is possible to produce models lacking the leading "OG" pulse. This picture most easily produces pulses with ∼180° phase separations and tends to have "Off-pulse" flux comparable to the "Bridge" between the two peaks. In its original form this Two-Pole Caustic (TPC) model (Dyks & Rudak 2003) truncates emission at the lesser of r_LC from the star or a distance of 0.75r_LC from the spin axis. This basic model uses the same w relation as the OG model (Equation (5)).

Modifications of this geometry have been suggested as better approximations of the physical realization of such vacuum zones in the Slot Gap model (Muslimov & Harding 2004; Grenier et al. 2010). For the first version (TPC2), the emission extends to 0.95r_LC and utilizes a slower, plausibly more physical, w scaling:

$$\begin{equation} w=\case{1}{2}\,({10^{33}\,{\rm erg\ s}^{-1}}/{\dot{E}})^{1/6}, \end{equation} \tag{ 6 }$$

which has wide w > 0.1 gaps even for very energetic Crab-like pulsars. Finally, we also define a full gap version of this model (TPC2_FG), with γ-ray emission coming from the full volume between the last closed field lines and the w value calculated from Equation (6).

There are also now a promising set of non-vacuum models, generated by numerical realizations of the force-free magnetosphere. One version, the "Separatrix Layer" (SL) model (Bai & Spitkovsky 2010) generates emission from field lines that start from w ≈ 0.1–0.15 and merge in the wind zone. In this picture, γ-ray emission arises from near the light cylinder, extending several tens of percent past this distance. While this model has some attractive features, especially for young pulsars with robust pair production, it requires numerical realizations. We concentrate here on comparing the Fermi sample with the analytically derived vacuum models (OG and TPC) and defer comparison with the SL model to future work.

The most uniform γ-ray pulsar sample available at present is the set of 38 young objects reported in the First Fermi LAT Catalog of γ-Ray Pulsars (Abdo et al. 2010a). We will refer to this as the "6 Month" sample (because that work was based on six months of Fermi LAT data). We also consider a larger sample of 50 young γ-ray pulsars, including all non-recycled pulsar detections announced as of 2010 July. In particular, this includes the "blind search" detections described in Saz Parkinson et al. (2010), as well as several individually announced discoveries. We refer to this as the "Current" sample; typically ∼1 year of LAT exposure contributed toward these detections.

Fermi has two different routes to γ-ray pulsar discovery. When the signal is folded on an ephemeris, based typically on radio observation, significant pulsations can be seen to quite low levels. We term these "radio-selected" pulsars. Pulsations can also be discovered "blindly," by searching directly in the GeV photons, but this method requires a somewhat higher flux for detection. We use here a threshold increase of 3× for such "γ-selected" detections (Abdo et al. 2010a).

A simple division of pulsars into these two classes is complicated by the discovery of several γ-ray pulsars with radio luminosities an order of magnitude or more fainter than those of the typical pulsar population (Camilo et al. 2009; Abdo et al. 2010b). These sources were detected in deep, targeted integrations of young supernova remnants, unidentified γ-ray sources or detected γ-ray pulsars. Their low flux density would not be accessible in the typical sensitivity of modern large area sky surveys, 0.1 mJy at 1.4 GHz.

Since we are simulating the properties of sources detectable in large sky surveys, we wish to assign real detected pulsars to one of these classes based on their observed fluxes, not the accident of whether they had exceptionally deep radio observation. Thus, "radio-selected" objects must have the survey-accessible flux density above. γ-selected objects must have the larger GeV flux required for "Blind search" discovery. This results in a re-classification of several objects. There are three pulsars that had prior radio detections with 1.4 GHz flux densities below 0.1 mJy (PSR J0205+6449, 0.04 mJy; PSR J1124−5916, 0.08 mJy; and PSR J1833−1034, 0.07 mJy). For the purposes of this work, we do not consider these objects to be radio-selected γ-ray pulsars. Two of the three objects have high enough γ-ray flux levels that even without the assistance of a radio ephemeris-guided folding search, they would still have been detected as γ-selected pulsars (PSR J0205+6449 and PSR J1833−1034). Thus, these two pulsars are added to the γ-selected subset of the observed population. The third object, PSR J1124−5916, likely could not have been detected to date in blind γ-ray searches, and so is removed from the observed population. Finally, PSR J2032+4127 was found in a blind γ-ray search, but radio follow-up observations found a 1.4 GHz flux of 0.24 mJy. Thus, this object "should have" been detected in radio surveys and is added to the radio-selected subset of the observed population.

With these classification amendments, the "6 Month" sample of 37 young γ-ray pulsars contains 19 radio-selected and 18 γ-selected pulsars. Similarly, the "Current" sample contains 24 radio-selected and 26 γ-selected objects. These are the samples against which we will test the γ-ray emission models.

In our simulation, any γ-ray pulsar with a radio beam crossing the Earth line of sight with modeled flux density >0.1 mJy is assumed detectable in a survey and termed radio-selected. For the γ-ray pulsars, we model the phase averaged flux on the Earth line of sight and compare with the estimated detection threshold for each sky location plotted as a map in Abdo et al. (2010a). We scale the ephemeris detection threshold with the square root of exposure time from the six-month estimates in that paper. To be γ-selected, a modeled pulsar must exceed 3× this ephemeris detection threshold, and must have a modeled radio flux density <0.1 mJy at Earth.

The first characteristic typically measured for a light curve is the number of peaks in the pulse profile. To compare with the data, we have taken each pulsar in our simulated Galaxy, applied the beam shape for each of the five trial γ-ray models described above, and tagged peaks with an algorithm that identifies peaks and their phases in the modeled γ-ray light curves. To increase similarity to the actual data we blocked the model light curves into 50 phase bins, before running the peak tagging algorithms. The results range from 0 peaks (γ-ray flux on the Earth line of sight, but no strong peak detected) to 4 or more distinct sharp peaks. Since for some models the peak morphology varies significantly with pulsar age/spin-down luminosity, we divide the sample into four bins of ${\dot{E}}$, logarithmically spaced.

The results are shown as histograms in Figure 4, with bars for each of the five models. For comparison, the measured multiplicities (always 1 or 2) for the "Current" sample are shown by the circles in each panel. It is immediately apparent that the original TPC model (central bars) produces many three and four peak light curves not seen in the data. The amendments to a slower w evolution (Equation (6)) and gaps extending to 0.95R_LC (wide right bars) make this disagreement worse. Extending the emission throughout the full gap (narrow right bars) does blur out some weaker peaks, mitigating the problem. We will use this blurred version (slow w evolution, fully illuminated gaps) in the remainder of the paper when we refer to "TPC," as it and the original TPC model give the best match to the observed quantities. The differences in the OG predictions between the Gaussian illumination (wide left bars) and the full gap illumination (narrow left bars) is relatively minor.

Zoom In Zoom Out Reset image size

We quantify the data comparison in Tables 1 and 2. Rows 1 through 4 of Table 1 show the Poisson probability that the model is a fit to the Fermi data. Each row represents one energy bin, and the best-fit model for that bin has its probability in bold. This emphasizes that the TPC2 model fairs particularly poorly without the full gap illuminated, and that for each bin the OG models provide substantially better representations of the peak multiplicity. The modeled presence of 0 peak sources (which could not, by definition, be discerned in the data) lowers the probability for all models. Also, the absolute normalization is meaningful here, as these computations are run with the same model-determined pulsar birthrate (see Section 6). In row 5, we calculate the probability decrement of each model when compared to the best model (original OG), summed across the four energy bins. Again, in the rest of the paper we adopt the original OG and TPC2_FG version as the representatives of the two fiducial models.

Table 1. Model Probabilities from Peak Number Comparison

$\dot{E}$	log(Poisson Probabilities)
	OG	OGFG	TPC	TPC2	TPC2FG
${\rm LOG}(\dot{\rm E})\ > \ 36.5$	−4.46	−4.66	−5.58	−8.42	−6.02
$36.5\ > \ {\rm LOG}(\dot{\rm E}) \ > \ 35.5$	−3.17	−3.59	−5.82	−11.22	−7.38
$35.5\ > \ {\rm LOG}(\dot{\rm E}) \ > \ 34.5$	−2.54	−2.75	−5.52	−8.79	−4.95
$34.5\ > \ {\rm LOG}(\dot{\rm E}) \ > \ 33.5$	−3.43	−3.19	−4.97	−7.03	−5.80
ΣΔLog(P)	...	−0.59	−8.29	−21.86	−10.55

Download table as:  ASCIITypeset image

Table 2. Model Probabilities from Pulse Morphology Comparisons

$\dot{E}$	Δ–δ : 2D K-S		Δ : 1D K-S		Bridge/Off-pulse : FoM
	OG	TPC2FG	OG	TPC2FG	OG	TPC2FG
${\rm LOG}(\dot{\rm E})\ > \ 36.5$	−0.96	−3.40	−2.00	−2.05	600 ± 128	462 ± 142
$36.5\ > \ {\rm LOG}(\dot{\rm E}) \ > \ 35.5$	−0.92	−3.52	-0.97	−0.72	1928 ± 268	387 ± 143
$35.5\ > \ {\rm LOG}(\dot{\rm E}) \ > \ 34.5$	−2.30	−4.52	-1.18	−0.83	552 ± 198	316 ± 159
$34.5\ > \ {\rm LOG}(\dot{\rm E}) \ > \ 33.5$	...	...	−0.65	−1.52	221 ± 138	143 ± 116
ΣΔLog(P)	...	−7.26	...	−0.32	...	−7.37

Download table as:  ASCIITypeset image

A particularly powerful test of the models arises from comparing the γ-ray peak separations Δ (available for all light curves with two or more peaks, taken as the widest strong peak separation) and the lag of the first γ-ray peak from the magnetic dipole axis δ. For radio-detected pulsars this axis is typically marked by the main radio pulse. The "Δ − δ" plot of these quantities is an excellent probe of the model geometry (Romani & Yadigaroglu 1995; Watters et al. 2009). We produce such plots here in Figure 5, again with four different $\dot{E}$ bins.

Zoom In Zoom Out Reset image size

γ-ray pulsars that are not detected in the radio will have an undefined value of δ; accordingly these objects are plotted in a band to the left, beyond the vertical dashed line at δ ∼ −0.06. Similarly, pulsars with only a single peak in their γ-ray pulse profiles have an undefined value of Δ. We assign to these singly peaked profiles a value Δ = 0 and plot in a band along the x axis. Note that in practice limited signal-to-noise ratio (S/N) in the γ-ray profiles usually prevents measurements of peak separations any closer than Δ ≈ 0.2, so the observed (large dots) single pulse profiles doubtless include some unresolved tight doubles, such as appear in the models particularly for small ${\dot{E}}$.

In the plots, the data show a clear inverse correlation between δ and Δ. This was predicted for outer magnetosphere models (Romani & Yadigaroglu 1995) and indeed appears for the (dark) OG model points. The trend is weak or absent in the (light) TPC points. In the OG case the correlation is easily understood since both caustic peaks are generated from the magnetic pole opposite to that producing the radio emission. As the Earth line-of-sight cuts a smaller chord across this hollow cone, δ increases and Δ decreases. In contrast, the TPC model produces the first γ-ray peak from the same pole as the radio emission. The δ is nearly fixed, showing only small variation with viewing angle and age.

A less obvious difference is the OG correlation of Δ and $\dot{E}$, again weak or absent in the TPC model, which as expected shows a typical Δ = 0.4–0.6, with a strong concentration at Δ ≈ 0.5, for all spin-down luminosities. In the LAT data, if we consider for each energy bin the fraction of objects with Δ>0.4 we find 73% (11/15), 56% (9/16), 38% (5/13), and 17% (1/6), moving through our four energy bins from high to low.

For a quantitative measure of the goodness of fit between the data and the models, we employ a two-dimensional Kolmorgorov–Smirnov test (K-S test) in this plane on the radio-selected samples. We can also run a standard one-dimensional K-S test on the γ-selected objects, for which only a Δ value exists. The probabilities from these tests are displayed in Table 2. With only a single detection, the two-dimensional K-S test is unavailable for the lowest $\dot{E}$ bin. The bottom row of Table 2 (as for Table 1) gives the summed probability decrement compared to the best (OG) model; although the OG model itself is not adequate (Prob ∼10⁻⁴), probabilities for the TPC models are factors of ⩾10⁷ lower. Statistically, the OG model is strongly preferred.

Single peaked (x-axis) γ-ray light curves deserve some additional discussion. As noted above, objects may appear here when the data do not resolve a double. Such objects in the OG picture are expected to appear at δ ≈ 0.3–0.4. In the OG model, a peak may also be single if the first γ-ray peak is missing. Since this caustic forms at high altitudes, it is rather sensitive to field line distortion from sweep back and to aberration effects, especially for moderate to large $\dot{E}$. In the vacuum models, we find that this caustic (and peak) are often missing for energetic pulsars when the observer line of sight lies well away from the spin equator. Near the spin equator the caustics are unaffected, so that large Δ appear, but smaller Δ < 0.3 are missing for energetic pulsars. As a result δ is large (∼0.5–0.6) since it now measures the distance to what is normally the second pulse. This effect may be seen in the models and data of the left panels of Figure 5. However, it should be noted that the sensitivity of the first peak caustic to field line perturbations makes the extent of this effect difficult to predict. For example, Romani & Watters (2010) found that open zone currents can reduce or eliminate such first peak "blowout," so realistic models containing plasma may affect the prevalence of large δ.

Finally, it should be noted that if small altitude emission occurs, then the trailing side of the radio-producing pole may form a caustic (this is the first peak in the TPC picture). As seen by the heavy concentration of TPC model dots in Figure 5, this results in δ ≈ 0.1. The Fermi LAT may have in fact uncovered one such object: PSR B0656+14, which appears near δ = 0.2 in the third panel, has a single pulse, a peculiar soft spectrum, and an apparent luminosity ∼30× smaller than seen from similar $\dot{E}$ pulsars. Polarization modeling suggests that it has a very small α and ζ, so that its outer magnetosphere beam should miss the Earth. In this interpretation, we only see the fainter low-altitude emission because of this pulsar's very low distance.

Thus, δ provides a useful way of sorting pulsar properties, even when only a single γ-ray peak is available. Improved S/N, polarization studies, and above all more physical magnetosphere modeling should help us calibrate this diagnostic.

The most common γ-ray profiles contain two peaks, and for these profiles we can introduce a third measure of pulse shape that can be applied to the bulk sample: the strength of the intra-peak ("bridge") and inter-peak ("off-pulse") emission. To estimate these flux levels we divide the double-peaked profiles into four phase windows. Two intervals of Δϕ = 0.14 (7 bins in a 50 bin light curve) are centered on the main peaks and the two remaining phase intervals (totaling 0.72 of pulse phase) are assigned to the "bridge" and "off-pulse." We next measure the mean flux in each of these phase windows. We then report the bridge fraction as 〈bridge〉/〈(bridge+peaks)〉 and the off-pulse fraction as 〈off-pulse〉/〈(bridge+peaks)〉. This easily defined measure generally corresponds well to a visual definition of bridge and off-pulse intervals, although it does not always isolate well the absolute pulse minimum. We plot these two fluxes for the model light curves, in the usual four panels of spin-down luminosity, in Figure 6.

Zoom In Zoom Out Reset image size

Clearly these quantities provide good model discrimination. Models radiating only in the outer magnetosphere (e.g., OG) have little off-pulse flux, especially for older pulsars, while models radiating at all altitudes have comparable bridge and off-pulse emission. In particular, the TPC picture produces essentially no light curves with off-pulse fluxes less than 40% of the pulse phase emission.

Unfortunately, unlike the peak number and separation, such fluxes are not provided in Abdo et al. (2010a), Saz Parkinson et al. (2010), and the other discovery papers. One difficulty in making such measurements is the identification of the true background level. In practice it can be very difficult to distinguish unpulsed "DC" magnetospheric emission from the contribution of an unresolved background in close proximity to the pulsar, such as expected from a surrounding pulsar wind nebula. Detailed study of source extension and spectrum at pulse minimum can help distinguish these cases, but this goes beyond the pulsar catalog data.

Nevertheless, the published light curves do plot an estimated background flux level. Accordingly we entered these light curves, defined the peak and intervening intervals exactly as above and measure the bridge and off-pulse flux levels. We assume a systematic 10% uncertainty in the reported background level, and then propagate this and the Poisson uncertainties in the bridge-, off-, and average-pulse fluxes to produce the two flux ratios and errors defined above. In a few cases the background level defined in Abdo et al. (2010a) was above the measured pulse minimum. In those cases, we reassigned the background flux to this level (and by definition the off-pulse flux was then 0). In a number of cases, the very large background pedestal prevented accurate measurements. While we show error bars for all estimates of double pulse pulsars in Figure 6, we mark with large dots only those measurements which had either a 2σ significance measure of both flux ratios or a ratio less than 0.1 at 2σ significance.

Bridge emission varies from ∼0.3× to >1× the average pulse flux at all spin-down luminosities, although a few high ${\dot{E}}$ objects seem to have virtually no emission away from the peaks. In contrast, the majority of the objects have 0.3× or less of the pulse emission in the off-peak window. Most of the best measured pulsars show very small values. Values of ∼0.1–0.3× appear for a number of the most energetic pulsars, but these measurements are of poor statistical significance. In any case such objects are most likely to show unpulsed contamination from associated PWNe, supernova remnants, and other diffuse, but unresolved emission. Thus overall, the general distribution of measured values matches better to the OG predictions, as shown in Table 2.

Nevertheless, several objects are distinctly inconsistent with emission from only high altitudes. In particular PSR J2021+4026 (panel 3) and PSR J1836+5925 (panel 4), show very strong off-pulse emission, which is spectrally consistent with magnetospheric emission (hard with a few GeV exponential cutoff). These unusual objects lie squarely within the TPC model zone. Interestingly both are quite low ${\dot{E}}$, and both have Δ very close to 0.5, consistent with a low-altitude, two-pole interpretation. Also, the very energetic PSR J1420−6048 (panel 1) shows strong off-pulse emission. Here, with Δ = 0.26 a low-altitude interpretation is not natural. However, unlike the other two cases, the very large and poorly determined background (from the bright PWN emission in the surrounding "Kookaburra" complex) may exceed the catalog background flux level. However, it is clear that, at least in a few cases, an emission component other than the outer magnetosphere of the OG model contributes substantially to the pulse profile. An excellent candidate is low-altitude emission, such as posited in the TPC model. An alternative is a current-induced shift of the OG start below r_NC (Hirotani 2006). It will be particularly interesting to discover the physical effect that causes such emission to be significant for a few pulsars, but negligible for most.