Searching a Thousand Radio Pulsars for Gamma-Ray Emission

The lowest photon fluxes among the 117 gamma-ray pulsars reported in the Second Fermi Large Area Telescope (LAT) Catalog of Gamma-ray Pulsars (Abdo et al. 2013, hereafter "2PC") were of order ${F}_{100}^{\min }\simeq {10}^{-9}$ photons cm⁻² s⁻¹, integrated above 100 MeV. The LAT's effective area averaged over the energy range where pulsars are bright is ∼5000 cm² (Atwood et al. 2009). Fermi's continuous all-sky survey covers every point on the sky with a roughly one-sixth duty cycle. Combining, we find that ${F}_{100}^{\min }$ yields a modest average of two photons per month recorded by the LAT, with >100% Poisson fluctuations in the monthly counting rate.

Detecting lower-flux pulsars gives access to a larger space volume, containing a rich variety of pulsars. Generally, luminosity scales with spindown power, ${L}_{\gamma }\propto \sqrt{\dot{E}}$ (see, for example, 2PC Figure 9). However, the dispersion around this trend spans over two orders of magnitude, presumably due to how the size, shape, and density of the emitting region depend on the neutron star's spin and magnetic field configuration. Luminosity may thus be intrinsically low for pulsars with rare combinations of properties, or may appear low because only a shoulder of the gamma-ray beam sweeps across the Earth (Romani et al. 2011). "Variety" therefore includes pulsars probing the $\dot{E}=4{\pi }^{2}I\dot{P}{P}^{-3}$ deathline for gamma-ray emission. P is the period, $\dot{P}=\displaystyle \frac{{dP}}{{dt}}$, and we set the moment of inertia to I = 10⁴⁵ g cm².

Pulsars may be "faint," i.e., near the LAT's detection threshold, even for high count rates: background is intense for low Galactic latitudes and in confused regions. Broad pulse widths also degrade sensitivity (Hou et al. 2014). Detecting and characterizing faint pulsars helps ensure that emission models are confronted with an accurate snapshot of the true pulsar population. It also helps determine the pulsar luminosity function needed to quantify the contribution of unresolved pulsars to the diffuse emission, for example, in the inner galaxy (Di Mauro et al. 2017).

We have been phase-folding LAT gamma-ray photons using hundreds of rotation ephemerides provided by the astronomers of the Pulsar Timing Consortium (Smith et al. 2008) since Fermi began routine observations on 2008 August 4 (MJD 54682). The Consortium provided ephemerides for ∼95% of known pulsars with spindown power $\dot{E}\gt {10}^{34}$ erg s⁻¹, and about half of those known with lower $\dot{E}$. This led to the discovery of gamma pulsations from roughly half of LAT's over 240 gamma-ray pulsars.¹⁴ The rest came from searches at the positions of pulsar-like unidentified LAT sources: a quarter of the total are radio-quiet pulsars found in blind gamma-ray periodicity searches (see Clark et al. 2017, and references therein), and nearly 90 are millisecond pulsars (MSPs) discovered via deep radio exposures (see, e.g., Pleunis et al. 2017). But flagging a source as pulsar-like means it is bright enough in gamma-rays to find its spectral shape and/or variability. The faint gamma-ray pulsars were all found by phase-folding with an ephemeris.

Weighting (Kerr 2011) is a key tool for these searches: the spatial and spectral map of the gamma-rays around the pulsar's position, combined with the LAT's energy-dependent point-spread function (PSF), translates into the probability of a given photon being signal or background (see Section 2). What is new in the present work is that we applied the simplified weighting method of Bruel (2018) to discover even fainter objects, for which a phase-integrated source is often not even detectable.

In Section 2 we summarize the Bruel (2018) photon weighting method. The new photon weighting scheme is fast, making it useful for brighter sources as well. In Section 3 we use the method to phase-fold 1269 pulsars, over 240 of which were known to pulse in gamma-rays as this work was being completed. The large sample highlights the stability and robustness of our search method, in particular, that the significance "noise floor" is safely below 4σ, lower than the 5σ threshold used in previous work. Gamma-ray pulsations from 16 previously known energetic radio pulsars are presented in Section 3. Folding the large sample also pushes the apparent gamma-ray emission "deathline" for pulsars to spindown power $\dot{E}\lesssim 8\times {10}^{32}$ erg s⁻¹ (Section 4). We summarize in Section 5.

A gamma-ray crosses the vacuum of space for thousands of years without incident, to enter perhaps into the LAT. Recoil by the atoms therein ensures conservation of momentum, allowing pair production, $\gamma \to {e}^{+}{e}^{-}$. The electron and positron ionize the silicon strips along their path in the LAT's tracker, but undergo multiple Coulomb scattering, mainly in the tungsten foils. Coulomb scattering causes the fundamental limit on the angular resolution of the LAT at low energy, and determines its PSF. The spacing between the silicon strips limits the minimum PSF above 10 GeV. Smith & Thompson (2009) provide further discussion.

At the end of the complex chain of LAT event reconstruction and selection, the LAT's energy-dependent PSF is modeled as a Moffat distribution (see Equation (3), below) in Δθ, the angular distance between a target direction (here, the pulsar's radio timing position) and a photon's sky direction $\vec{\theta }$. The distribution's width in degrees (68% containment) is

$$\begin{eqnarray}&&{\sigma }_{\mathrm{psf}}({E}_{\gamma })=\sqrt{{p}_{0}^{2}{\left({E}_{\gamma }/100\right)}^{-2{p}_{1}}+{p}_{2}^{2}},\end{eqnarray} \tag{ 1 }$$

where E_γ is the gamma-ray photon energy in MeV. For the P305 Pass 8 data¹⁵ used in this work (Atwood et al. 2013), the parameter values are p₀ = 5.445, p₁ = 0.848, and p₂ = 0.084.

To search for a pulsed gamma-ray signal, we select gamma-ray photons within 5^◦ of the pulsar position, approaching σ_psf for our E_γ > 50 MeV data sample. Their arrival times in the LAT are converted to a rotational phase of the neutron star using the parameters in the rotation ephemeris for that pulsar and the Fermi plug-in (Ray et al. 2011) to Tempo2 (Hobbs et al. 2006). We evaluate pulsation significance for the photon list with the H-test (de Jager & Büsching 2010). We stack (fold) the phases in histograms to examine and analyze the profile.

In the presence of background photons from nearby sources and diffuse emission, the signal-to-noise ratio can be optimized using different methods. The most effective to date is weighting (Kerr 2011). The weight $w({E}_{\gamma },\theta )$ is the probability that a photon with energy E_γ and apparent sky direction $\vec{\theta }$ comes from the pulsar rather than background. Kerr calculates the weights as $w={R}_{p}/({R}_{b}+{R}_{p})$, where R_p is the rate of gamma-rays from the pulsar, and R_b is the rate for gamma-rays from all the background components. The Fermi Science Tool gtsrcprob (see footnote 15) implements the method. It takes as input a "source model" describing the spectra and positions of all nearby gamma-ray sources, created using the gtlike tool, as well as ${\sigma }_{\mathrm{psf}}({E}_{\gamma })$ and $\vec{\theta }$, to calculate the expected rates and thus w. Kerr modified the H-test to include weights, summarized succintly in Section 3 of 2PC. Weighting for the 117 pulsars in 2PC, and for nearly all LAT team pulsar studies since 2011, was done with gtsrcprob.

However, the gtsrcprob tool has two significant drawbacks. First, gtlike requires high detection significance to allow gtsrcprob weights to be meaningful. The analysis can be complex for faint sources in highly confused regions, and may fail for faint sources where pulsations are nevertheless clearly detectable. Second, gtlike, and the tools gtdiffrsp and gtexposure that its use requires, impose much larger data files and longer computation times than needed for the simple weighting method.

2.1. Simplified Weighting

Bruel (2018) provides a simple formula that approximates the weights remarkably well,¹⁶

$$\begin{eqnarray}&&w({E}_{\gamma },{\rm{\Delta }}\theta )=g({E}_{\gamma },{\rm{\Delta }}\theta )\exp (-2{\mathrm{log}}^{2}({E}_{\gamma }/{E}_{\mathrm{ref}})),\end{eqnarray} \tag{ 2 }$$

with one free parameter for the energy scale, ${\mu }_{w}\,={\mathrm{log}}_{10}({E}_{\mathrm{ref}}/1\,\mathrm{MeV})$, discussed below. The geometrical factor g uses a Moffat distribution (Moffat 1969)¹⁷ to describe the angular distribution of the gamma-ray photons emanating from a point source,

$$\begin{eqnarray}&&g({E}_{\gamma },{\rm{\Delta }}\theta )={\left(1+\displaystyle \frac{9{\rm{\Delta }}{\theta }^{2}}{4{\sigma }_{\mathrm{psf}}^{2}({E}_{\gamma })}\right)}^{-2}.\end{eqnarray} \tag{ 3 }$$

The Moffat distribution resembles a Gaussian but with broader tails.

The exponential factor was obtained by studying weight distributions computed under the simple assumptions that the pulsar rate is much lower than the background rate, and that the background is uniform. Bruel (2018) explores the validity of these assumptions. Pulsar spectra are of the form

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dE}}_{\gamma }}=K{\left(\displaystyle \frac{{E}_{\gamma }}{1\mathrm{GeV}}\right)}^{-{\rm{\Gamma }}}\exp \left(-\displaystyle \frac{{E}_{\gamma }}{{E}_{\mathrm{cut}}}\right),\end{eqnarray} \tag{ 4 }$$

that is, hard spectra (typically 0.5 < Γ < 2.0) with sharp cutoffs (0.6 < E_cut < 6 GeV, see 2PC). K determines the flux. In the energy range germane to pulsars, the spectrum of the diffuse background is a steeper power law. In these conditions, the pulsar weights versus $\mathrm{log}{E}_{\gamma }$ are well-approximated by Gaussians with width 0.5, relatively independent of the pulsar's spectrum. What changes from one pulsar to another is the reference energy ${E}_{\mathrm{ref}}={10}^{{\mu }_{w}}$ for which the pulsar's weight distribution peaks. For a faint pulsar with a soft spectrum far off the Galactic plane (hence, with low background), E_ref can be below 1 GeV (μ_w = 3). A pulsar with a hard spectrum in a highly confused region can have a distribution peaking above 10 GeV (μ_w = 4).

2.2. Choosing μ_w

2.3. Rotation Ephemerides

In practice, we define H_max as the maximum H-test value for six trials: the three values of μ_w = 3.2, 3.6, and 4.0 described above, applied to two data samples. One sample uses 9.6 yr of LAT data, from the beginning of the Fermi mission, MJD 54682, until MJD 58182. The other uses LAT data only during the ephemeris validity period. Many pulsars rise to higher significance followed by a rollover at intermediate dates, but restricting to two epoch durations keeps the number of trials small and well-defined.

Figure 1 shows H_max versus $\dot{E}$ for 1269 pulsars, and Figure 2 shows a histogram of the same H_max values. A "noise floor" appears. We choose H_max > 25 (single trial weighted H-test significance of 4.1σ) as our threshold for detection of gamma-ray pulsations. With this threshold, we expect less than one false detection for our 1269 pulsar sample. We discuss the 16 gamma-ray pulsars thus discovered below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 also shows the distribution of H-test values obtained for the same pulsar sample, where the photon phases were replaced with uniform random values between 0 and 1, using μ_w = 3.6. We did this three times (solid histogram) and eight times (dotted histogram) for each pulsar and kept the maximum value, as in the procedure for H_max. The three-trial histogram matches that of the undetected pulsars. Three trials suffice because the μ_w = 2.3 trials are not included in Figures 1 and 2, but also because counting both data sets is redundant: the "all data" and "ephemeris validity" H-test values are in fact highly correlated. Two-thirds (one-third) of our ephemerides cover over half (over 90%) of the mission. Ephemeris validity is <500 days for only 5% of the sample. Eight trials adds the fourth scan, with μ_w = 2.3, for the two data samples, ignoring correlations between the samples. H_max > 25 never occurs in the random histograms.

A simple Monte Carlo calculation, illustrated in Figure 3, gives insight. Four simulated examples of significance growth versus time are shown. For the few photons per month signal rates evoked in the Introduction for pulsars near detection threshold, Poisson fluctuations can cause irregular growth, sometimes wildly so. We see such fluctuations in the data, and speculated that they could be due to changes in flux, instrument exposure, or ephemeris instability. In fact, Poisson fluctuations dominate. For the simulated signal and background rates in Figure 3, ${F}_{100}^{\min }\simeq {10}^{-9}$ photons cm⁻² s⁻¹, we simulated the Fermi mission 1000 times. The mean single-trial significance is ∼4σ, meaning that the pulsar is not detected half the time. The full width at half maximum of the significance distribution extends from 2.5σ to 5σ. We also used the simulation to understand how many years it takes to detect a pulsar having a given flux, duty cycle, and background environment. Statistical fluctuations smear the answer. Near threshold, the fluctuations are Poissonian and can be so large that H_max > 25 can occur a full year earlier or later than the average duration.

Zoom In Zoom Out Reset image size

3.1. Discoveries of Gamma Rays from Radio Pulsars

Table 1 lists 16 gamma-ray pulsars discovered using the above method. Twelve are young to middle-aged, and four are MSPs. Figures 4–6 show their phase-aligned >50 MeV gamma-ray and radio profiles, as well as the evolution of their H-test (pulsed Test Significance) since the start of the Fermi mission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Twelve Young (Top) and Four Recycled (Bottom) Radio Pulsars for which We See Gamma-Ray Pulsations

PSR	Ref.a	l	b	P	Distance	$\dot{E}/{10}^{34}$	S1400	H-test	Timing
		(°)	(°)	(ms)	(kpc)	(erg s−1)	(mJy)	signif.
J0729−1836	mlt+78	233.76	−0.34	510.2	2.4 / 2.9	0.56	1.4	41	P (J)
J1731−4744†	lvw68	342.56	−7.67	829.8	5.5 / 2.8 (0.7)	1.13	12.0	34	P
J1740+1000†	mca00	34.01	20.27	154.1	1.2 / 1.2	23.16	9.2	31	J (N)
J1757−2421	kom74	5.28	0.05	234.1	3.1 / 4.4	4.00	3.9	29	P (J, N)
J1816−0755	lfl+06	21.87	4.09	217.6	3.1 / 2.8	2.48	0.17	54	J
J1841−0524	hfs+04	27.02	−0.33	445.7	4.1 / 5.3	10.42	0.2	27	J (P)
J1853−0004†	hfs+04	33.09	−0.47	101.4	5.3 / 7.1	21.09	0.87	43	P (J, N)
J1913+1011	mhl+02	44.48	−0.17	35.9	4.6 / 4.8	287.14	0.5	32	J
J1925+1720	lbh+15	52.18	0.59	75.7	5.1 / 6.9	95.42	0.07	37	J
J1928+1746	cfl+06	52.93	0.11	68.7	4.3 / 5.8	160.39	0.28	31	J (N)
J1932+2220†	ht75b	57.36	1.55	144.5	10.9/ 7.5	75.38	1.2	44	J
J2208+4056†	slr+14	92.63	−12.09	637.0	1.0 / 0.8	0.080	0.48	50	J
J0636+5129	slr+14	163.91	18.64	2.86	0.2 / 0.5 (1.1+0.6−0.3)	0.569 ± 0.001		48	N
J1125−6014†	fsk+04	292.50	0.89	2.63	1.0 / 1.5	0.45 ± 0.3	0.05	35	P
J1327−0755	blr+13	318.38	53.85	2.68	25.0/ 1.7	0.36		25	N (G)
J1946+3417	bck+13	69.29	4.71	3.17	7.0 / 5.2 (2.5)	0.2 ± 0.1	0.29	29	J (E)

Note. The "Ref." column indicates the discovery paper. Distances are taken from the ATNF pulsar catalog. The first is derived from the dispersion measure (DM) using the YMW16 Model (Yao et al. 2017) except for J1932+2220, which uses the kinematic distance of Frail & Weisberg (1990), and for J0636+5129, for which Stovall et al. (2014) measured timing parallax. The second distance comes from the DM using NE2001 (Cordes & Lazio 2002). The additional distance for J0636+5129 comes from Arzoumanian et al. (2018) (note added in proof). The additional distance for J1731−4744 is from Shternin et al. (2017). The additional distance for J1946+3417 comes from the discussion of Figure 7, see the text. The rotation ephemeris used in this work (last column) was obtained from radio timing with N (Nançay Radio Telescope; Cognard et al. 2011); P (Parkes Radio Telescope; Weltevrede et al. 2010); J (Jodrell Bank Observatory; Hobbs et al. 2004b); E (Effelsberg; Barr et al. 2017); G (Green Bank Telescope; Kaplan et al. 2005). The observatories in parentheses also provided timing observations. The MSPs with $\dot{E}$ uncertainties are Shklovskii-corrected (see the text). The six pulsars with a † are colocated with LAT 8 yr sources.

^abck+13: Barr et al. (2013); blr+13: Boyles et al. (2013); cfl+06: Cordes et al. (2006); fsk+04: Faulkner et al. (2004); hfs+04: Hobbs et al. (2004a); ht75b: Hulse & Taylor (1975); kom74: Komesaroff, M. M. (1974, unpublished); lbh+15: Lazarus et al. (2015); lfl+06: Lorimer et al. (2006); lvw68: Large et al. (1968); mca00: McLaughlin et al. (2000); mhl+02: Morris et al. (2002); mlt+78: Manchester et al. (1978); slr+14: Stovall et al. (2014).

Download table as:  ASCIITypeset image

3.1.1. Nonrecycled Pulsars

The parameter space occupied by gamma-ray pulsars continues to fill in and expand. PSR J2208+4056 has $\dot{E}=8\times {10}^{32}$ erg s⁻¹, three times lower than the previous lowest spindown power for young gamma-ray pulsars. PSR J0729−1836 is now one of nine young pulsars with $\dot{E}\lt {10}^{34}$ erg s⁻¹. At the other end of the range, PSRs J1913+1011 and J1928+1746 both have $\dot{E}\gt {10}^{36}$ erg s⁻¹.

PSR J1731−4744 is the longest spin-period gamma-ray pulsar now known, with P = 830 ms. Slow pulsars (PSR J0729−1836 has P = 510 ms) attain $\dot{E}\simeq {10}^{34}$ erg s⁻¹ because their surface magnetic fields are high, B_S ≳ 10¹³ G. In contrast, PSR J1913+1011 is nearly as fast as the Crab but with B_S an order of magnitude smaller.

In the LAT 8 yr preliminary source list see footnote 15, FL8Y, only the six pulsars flagged in Table 1 lie within a source's 95% confidence ellipse. The rest are undetected, using a procedure similar to that of 3FGL (Acero et al. 2015), improved in preparation for 4FGL, the upcoming 4th Fermi LAT catalog. Weighted folding improves signal-to-noise roughly as the inverse of the pulse width δϕ, compared to a phase-integrated detection. Seven of the twelve young pulsars have Galactic latitude $| b| \lt 0\buildrel{\circ}\over{.} 6$ (and thus very high background levels), and all but one of these (PSR J1853−0004) have narrow gamma-ray pulses, with δϕ < 0.3 in phase.

The highest background levels are toward the Galactic center. The detection of PSR J1757−2421, only 53 away, is aided by its narrow pulse, δϕ < 0.1. PSR J1841−0524, with high $\dot{E}$ and intermediate distance, presumably needed over nine years of LAT data to be seen because its overlap with the bright pulsar wind nebula HESS J1841−055 raises its background (Acero et al. 2013). Only 09 separate high $\dot{E}$ PSRs J1925+1720 and J1928+1746, raising each's background marginally compared to the already high level of this busy region in the Galactic plane, compounding detection difficulty.

Finally, PSR J1932+2220 is the most distant of the new pulsars. The NE2001 (Cordes & Lazio 2002) and the YMW16 (Yao et al. 2017) Dispersion Measure (DM) models place it at 7.5 and 8 kpc from Earth, respectively, whereas its distance derived from H i absorption is ${10.9}_{-0.8}^{+0.3}$ kpc (Verbiest et al. 2012). The others pulsars lie in the range of ∼1 to ∼6 kpc, like the bulk of the 2PC pulsars. PSR J1932+2220's "faintness" is thus a combination of low flux (medium $\dot{E}$ at large distance) and the background level just off the plane.

3.1.2. MSPs

Three of the four new MSPs are "faint" simply due to low gamma-ray fluxes: neither the intense background at low latitudes nor broad peaks made them hard to detect. Three have distances d < 2 kpc, as derived from their DM values. The only one in the FL8Y list, PSR J1125−6014, is also the only one deep in the plane, with b = 09, with five other gamma-ray pulsars less than 2° away. It also has the narrowest pulse, facilitating detection. In Figure 6 its signal is contained in a single phase bin of width δϕ = 1/50 = 0.02 in phase, so its pulse is <53 μs wide.

Barr et al. (2017) used timing of the binary orbit and radio polarization measurements to characterize PSR J1946+3417. They argue that the more intense radio peak (see Figure 6) could come from a polar cap, while the smaller one would come from the outer magnetosphere. They determined the pulsar's geometry to be favorable for gamma-ray viewing from Earth, and attributed its non-detection by the LAT to a large distance and low spindown power. Their reasoning overlaps that invoked for young radio-interpulse pulsars in Section 4.1. PSR J1946+3417 is shown by the lowest large triangle in Figure 7, discussed below.

Zoom In Zoom Out Reset image size

Having now detected it in gamma-rays, we consider Barr et al's prediction that if indeed radio emission comes from an outer gap (OG), the gamma-ray emission would be aligned in phase with the second radio peak. The gamma-ray profile for PSR J1946+3417 in Figure 6 shows two peaks separated by 0.25 in phase. The radio peak they attribute to the OG is between them, on the inner shoulder of the trailing gamma peak, consistent with their prediction. The uncertainty shown is derived from the difference between the DM value in our ephemeris and that published by Barr et al. (2017). A wealth of accurate observational information now exists for this pulsar, making it an interesting candidate for detailed modeling.

For any pulsar, requiring gamma-ray efficiency $\eta \,={L}_{\gamma }/\dot{E}\lt 100 \%$ reveals its maximum possible distance, for a given assumption about its beaming fraction f_Ω, because the luminosity ${L}_{\gamma }=4\pi {d}^{2}{f}_{{\rm{\Omega }}}{G}_{100}$ depends on the distance and on the energy flux integrated above 100 MeV, G₁₀₀. For MSPs, the Shklovskii effect makes the η < 100% constraint even stronger, because the intrinsic (corrected) spindown power ${\dot{E}}^{\mathrm{int}}$ is also sensitive to distance, as illustrated in 2PC Figure 11. We now explore the plausible distances of the four MSPs, and hence their possible ${\dot{E}}^{\mathrm{int}}$ values.

The Doppler correction to $\dot{P}$ depends on the system's proper motion transverse to the line of sight, μ (Shklovskii 1970), and the acceleration due to its and the Sun's locations the Galaxy. From the observed spindown $\dot{P}$, the intrinsic value is ${\dot{P}}^{\mathrm{int}}=\dot{P}-{\dot{P}}^{\mathrm{shk}}-{\dot{P}}^{\mathrm{gal}}$, with

$$\begin{eqnarray}&&{\dot{P}}^{\mathrm{shk}}=\displaystyle \frac{1}{c}\,{\mu }^{2}\,d\,P=k{\left(\displaystyle \frac{\mu }{\mathrm{mas}{\mathrm{yr}}^{-1}}\right)}^{2}\left(\displaystyle \frac{d}{\mathrm{kpc}}\right)\left(\displaystyle \frac{P}{{\rm{s}}}\right)\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\dot{P}}^{\mathrm{gal}}=\displaystyle \frac{1}{c}\,{{\boldsymbol{n}}}_{10}\cdot ({{\boldsymbol{a}}}_{1}-{{\boldsymbol{a}}}_{0})P\end{eqnarray} \tag{ 6 }$$

where k = 2.43 × 10⁻²¹. The Galactic potential model of Carlberg & Innanen (1987) and Kuijken & Gilmore (1989) provides the accelerations ${\boldsymbol{a}}$₁ of the pulsar and ${\boldsymbol{a}}$₀ of the Sun. Here, ${\boldsymbol{n}}$₁₀ is the unit vector in the direction from the solar system to the pulsar.

For PSR J0636+5129, the two terms nearly cancel and the Shklovskii correction is small, <1%, as are the uncertainties, owing to a timing parallax distance measurement and a small, precise μ value (Stovall et al. 2014). PSR J1327−0755, in contrast, lies on a line of sight where converting the DM to distance is highly uncertain (see Table 1) and its nominally large μ value has large uncertainties. We therefore do not correct $\dot{P}$ for this MSP. The larger timing parallax distance found by Arzoumanian et al. (2018) does not change this result qualitatively.

For PSR J1125−6014 our ephemeris has μ = 17 ± 0.1 mas yr⁻¹. Of the two distances in Table 1, the closer one, d = 1 kpc, gives a correction of ${\dot{E}}^{\mathrm{int}}=0.6\dot{E}$. The 3FGL value of G₁₀₀ (Acero et al. 2015) and f_Ω = 1 then give a typical value of η = 30%. The farther distance implies small f_Ω to reach that range of η values.

Since PSR J1946+3417 is undetected in FL8Y, we have no G₁₀₀ value, and thus no η < 100% distance constraint. However, the "heuristic" energy flux ${G}_{h}=\sqrt{\dot{E}}/(4\pi {d}^{2})$ can also flag suspicious distances. In Figure 7, G_h versus $\dot{E}$ for our twelve hundred pulsars, MSP J1946+3417 is well below any other gamma-ray pulsar, implying that the two DM distances in Table 1 are overestimated. To raise PSR J1946+3417 to G_h = 8 × 10¹⁴ (erg ${{\rm{s}}}^{-1}{)}^{1/2}$ kpc⁻² requires d = 2.5 kpc, two or three times closer than the DM distances.

The Shklovskii correction lowers G_h for PSR J1946+3417 even more. Using μ = 8.64 ± 0.25 mas yr⁻¹ from Barr et al. (2017), for d = 7 kpc, observed $\dot{E}=3.8\times {10}^{33}$ erg s⁻¹ corrects to a record-breaking ${\dot{E}}^{\mathrm{int}}=0.83\times {10}^{33}$ erg s⁻¹. For $d=5.2\,\mathrm{kpc}$, the corrected value is as small as any seen for a gamma-ray MSP, ${\dot{E}}^{\mathrm{int}}=1.8\times {10}^{33}$ erg s⁻¹, while d = 2.5 kpc suggested by Figure 7 gives a more typical value of ${\dot{E}}^{\mathrm{int}}=2.8\times {10}^{33}$ erg s⁻¹. We add the d = 2.5 kpc to Table 1 to underline that the gamma-ray detection favors a closer distance. Finally, Barr et al. (2017) determined this pulsar to have 1.828 ± 0.022 solar masses, making I, and thus $\dot{E}$, 30% larger than obtained using the usual Chandrasekhar mass for the neutron star. This shifts how a gamma-ray flux constrains this pulsar's distance.

3.2. The Birthplace of PSR J1731−4744

Striking in Figure 1 is a minimum spindown power $\dot{E}$ for gamma-ray pulsars. Most with $\dot{E}\lt {10}^{34}$ erg s⁻¹ are MSPs. The discovery of gamma-ray pulsations from PSR J2208+4056 establishes that while it may be rare for such low-power young pulsars to be seen with the LAT, it is not impossible. We discuss this pulsar further below.

This new sample allows checks on attempts to understand which pulsars will or will not be seen in gamma-rays. Rookyard et al. (2017) focused on the radio profiles of young pulsars with $\dot{E}\gt {10}^{35}$ erg s⁻¹, with the idea that the profiles result from the inclinations of the pulsar's rotation and magnetic axes relative to the line of sight (ζ and β, respectively). We have seven such pulsars in Table 1, five of which are in their sample. In their Figure 2, radio pulse width w₁₀ versus $\dot{E}$, gamma-ray detected pulsars lie mostly below a diagonal ${w}_{10}=48\mathrm{log}\dot{E}-1674$, for $\dot{E}$ in units of erg s⁻¹. Our seven straddle this line. Nevertheless, the general trend they pointed out—that the third of high-$\dot{E}$ pulsars that remain undetected in gamma-rays tend to be those with broader radio pulses—is borne out.

For the MSPs, Guillemot et al. (2016) used timing parallax distances, with well-determined uncertainties, to confirm that once relative motion effects have been taken into account via the Shklovskii proper motion corrections, the observed $\dot{E}\approx {10}^{33}$ erg s⁻¹ MSP deathline value is robust. The four new gamma-ray MSPs reported here are consistent with this.

Again for MSPs, Guillemot & Tauris (2014) also look for observations that constrain the inclination angles, to then predict whether a gamma-ray beam might sweep the Earth. During the MSP recycling process (i.e., during spin-up), the mass of the precursor of the white dwarf companion determines the orbital radius, and thus the orbital period P_B. Using P_B and ${a}_{p}\sin i$, they deduce i, the orbital inclination. They use i ≈ ζ, and expect that gamma-rays can be seen mainly for large ζ, because emission models generally favor equatorial gamma-ray beams. Three of our four MSPs have P_B > 1 day, the limit for their model's validity. For two, the (P_B, a_p sin i) values lead to ζ values near 45°, and for the third ζ ≈ 70°. These values fill in the ζ histogram of their Figure 5, tending to support their conclusions.

4.1. PSR J2208+4056

4.2. Distance Biases

4.3. Spectra

In its eleventh year in orbit, Fermi/LAT continues to discover about 24 gamma-ray pulsars per year. Two discovery channels are deep radio searches and blind gamma-ray periodicity searches at the positions of bright pulsar-like unidentified gamma-ray sources. This paper concerns a third channel, phase-folding gamma-rays using a rotation ephemeris obtained from sustained observations of previously known pulsars. This channel is the most sensitive to gamma-faint pulsars.

Gamma-ray pulsations for 16 pulsars were presented, 12 young and 4 recycled, made possible because (i) several hundred updated (1269 in total) ephemerides were applied to a 9.6 yr gamma-ray data set; (ii) Bruel's (2018) gamma-ray weighting technique can be applied regardless of whether or not the pulsar is seen as a phase-integrated point source; (iii) we demonstrated that the detection threshold can be lowered to 4.1σ without incurring false positive detections. The new gamma-ray pulsars include the faintest yet seen, giving fewer than 10 photons per hundred days. We described some of the consequences of the resulting large Poisson fluctuations. PSR J2208+4056 has $\dot{E}=8\times {10}^{32}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, three times lower than what had been the lowest spindown powers for known gamma-ray pulsars these last several years.

We discussed selection biases associated with the new discoveries. Most of the new pulsars have very narrow peaks, improving sensitivity, and lie at low Galactic latitude, where the background is intense. The Third Gamma-ray Pulsar Catalog ("3PC," in preparation) will contain more than twice as many gamma-ray pulsars as 2PC. The unambiguous identification and localization of gamma-ray pulsars helps improve 4FGL, the upcoming fourth Fermi/LAT catalog of sources using 8 yr of data. Both catalogs should go to press in the coming year.

Emission modelers should bear in mind that because faint pulsars are easier to detect if the pulses are narrow, tuning models to match the observed sample would be a mistake. If, for example, pulses broaden as $\dot{E}$ decreases, those pulsars would be absent from the current observed sample.

We discussed whether the rarity of detected gamma-ray pulsars with spindown power $\dot{E}\lt {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ is a feature of the gamma-ray emission mechanism (a "deathline"), or if it is due to observational selection effects. There does seem to be a deathline, with luminosity dropping below the ${L}_{\gamma }\propto \sqrt{\dot{E}}$ trend at low $\dot{E}$. But faint emission at low $\dot{E}$ for rare pulsars may be undetected at present because the distances implied by the large volumes required to obtain a pulsar sample large enough to provide a few detections impose values of ${G}_{h}=\sqrt{\dot{E}}/(4\pi {d}^{2})$ below the LAT's sensitivity. Sensitivity is enhanced for narrow pulses and for high Galactic latitudes, but requiring those conditions decreases the sample size. We also mentioned model predictions that peak emission would shift from the GeV to the MeV energy range as $\dot{E}$ decreases, making LAT detections more difficult.

Along with the simple scheme applied here, Bruel (2018) presents a pulsar weighting method that is more sensitive but also more difficult to use. Future work will be to refold the thousand pulsars, as well as other pulsars for which ephemerides may become available, with the refined method. Preliminary indications are that we will find still more faint gamma-ray pulsars.

The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS).

The Parkes radio telescope is part of the Australia Telescope, which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO.

The Lovell Telescope is owned and operated by the University of Manchester as part of the Jodrell Bank Centre for Astrophysics with support from the Science and Technology Facilities Council of the United Kingdom.

The Robert C. Byrd Green Bank Telescope (GBT) is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

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 performed in part under DOE Contract DE- AC02-76SF00515. Work at NRL is supported by NASA.