DISCOVERY OF NINE GAMMA-RAY PULSARS IN FERMI LARGE AREA TELESCOPE DATA USING A NEW BLIND SEARCH METHOD

The Fermi Gamma-ray Space Telescope has been operating since its launch in 2008 June. The Large Area Telescope (LAT) on board the Fermi satellite has an effective area of ∼0.8 m² (on-axis, above 1 GeV) and is sensitive to photons with energies from 20 MeV to more than 300 GeV (Atwood et al. 2009). Post-processing assigns an arrival time, energy E, and direction to the photons; we call these "events." The arrival times are accurate to better than 1 μs, and the energies are accurate to within 15% between 0.1 and 10 GeV on-axis. The directional precision is energy dependent: 68% of photons have angular offset less than ∼08 × (E/GeV)^−0.8 from the true direction (Abdo et al. 2009c).

Gamma-ray pulsars are among the most interesting sources observed by the Fermi LAT. These are rapidly spinning neutron stars whose regular beam of gamma-ray emissions passes by the detector with each rotation. The Fermi LAT has detected gamma-ray pulsations from more than 50 normal and millisecond pulsars (MSPs) previously discovered in radio-frequency searches (see, e.g.; Abdo et al. 2010c; Ransom et al. 2011; Keith et al. 2011; Cognard et al. 2011; Theureau et al. 2011).

In contrast, the so-called blind searches for gamma-ray pulsars are not guided by any prior knowledge (e.g., from radio or X-ray observations) of the pulsars' parameters. Previous blind searches of the data recorded by the Fermi-LAT have been spectacularly successful; they have discovered 26 gamma-ray pulsars (Abdo et al. 2009a; Saz Parkinson et al. 2010; P. M. Saz Parkinson et al. 2011, in preparation). This paper describes the "blind" discovery of nine additional pulsars using a new method.

The blind search problem is computationally demanding because the relevant pulsar parameters (typically sky position, frequency f, and spin-down rate $\dot{f}$) are not known a priori and must be explicitly searched (Chandler et al. 2001). So far, blind searches of LAT data have used a clever "time-differencing technique" as described in Atwood et al. (2006) and Ziegler et al. (2008). One powerful motivation for seeking even better methods is the application to MSPs, of which previous blind searches have not found any. The blind search for MSPs is more computationally challenging due to the higher frequency range that must be covered. Moreover, most MSPs are in binary systems, where the orbital modulation parameters must also be searched, adding orders-of-magnitude to the computational complexity and challenge.

This paper presents results from a new effort to find isolated gamma-ray pulsars (including MSPs), using a novel method inspired by computationally efficient techniques recently developed to search gravitational-wave detector data for weak continuous-wave signals from rapidly spinning isolated neutron stars (Pletsch & Allen 2009; Pletsch 2010). In particular, the search method uses a recently developed optimal incoherent-combination method first described in Pletsch & Allen (2009) together with a "sliding coherence window technique" (Pletsch 2011).

The method was originally intended to find isolated MSPs up to 1.4 kHz spin frequency with spin-down rates in the range −5 × 10⁻¹³ Hz s$^{-1} \le \dot{f} \le 0$, with a corresponding range of characteristic ages $\tau =-f/2\dot{f}$. However, the search technique is also sensitive to normal (non-MSP) isolated pulsars, and several normal pulsars were discovered soon after we began. This prompted us to extend the spin-down-rate search range, in order to include younger objects (having smaller τ), down to ages (∼ kyr) comparable to the Crab pulsar. Here we report on the discovery of nine pulsars from this ongoing effort in which a total of 109 sources selected from the Fermi-LAT Second Source Catalog (Abdo et al. 2011b) are being searched for new gamma-ray pulsars.

The outline of this paper is as follows. Section 2 describes the LAT data preparation and the selection of unidentified sources to search for previously unknown gamma-ray pulsars. The new hierarchical search method is explained in Section 3 and illustrated in Section 4 with a detailed example: our first discovery, PSR J1620–4927. Section 5 presents the results for all the new gamma-ray pulsars. The search for counterparts in other regions of the electromagnetic spectrum is discussed in Section 6. This is followed by a brief conclusion.

The Fermi-LAT Second Source Catalog (2FGL; Abdo et al. 2011b) lists 1873 sources, described by fits to elliptically shaped 95% confidence sky regions. Among these sources 576 are not associated with counterparts observed at other wavelengths and thus might contain unknown gamma-ray pulsars.

In searching for new gamma-ray pulsars, it is important to identify and exclude sources that are blazars. Previously observed Fermi-LAT gamma-ray pulsars (see, e.g., Abdo et al. 2010c) have sharp cutoffs in their emission spectra at a few GeV, and stable gamma-ray fluxes. In contrast, blazars emit above 10 GeV, and their fluxes vary with time. An illustration of the different spectral properties and variability behavior can be seen in Figure 17 of Abdo et al. (2011b). Gamma-ray pulsars tend to have large curvature significances ("Signif_Curve" in the 2FGL catalog, which gives the improvement in the quality of the spectral fit when changing from a power law to a curved spectral model) and small variability indices ("Variability_Index" in the 2FGL catalog, a measure of flux instability over time). In contrast, blazars tend to have small curvature significances and large variability indices.

For the search for gamma-ray pulsations we select 2FGL sources with curvature-significance values greater than 4σ (here and throughout the manuscript, σ denotes the standard deviation for a Gaussian distribution) and variability indices smaller than 41.6. We further select bright objects by choosing sources with detection significances ("Signif_Avg") greater than 10σ. Finally, we restrict the list to sources with no known associations ("unassociated sources") or associated with known supernova remnants (SNRs). By applying the above selection criteria to a preliminary version of the 2FGL catalog, we obtained a list of 109 2FGL sources to search for gamma-ray pulsations. (One source, 2FGL J0621.9+3750, was associated with an active galactic nucleus (AGN) in the final 2FGL catalog.)

The Fermi Science Tools (ST)²⁴ v9r23p1 are employed to select the Fermi-LAT events for our search. Using gtselect, we take events from 2008 August 4 to 2011 April 6, with reconstructed directions within 8° of the gamma-ray sources, energies above 100 MeV, and zenith angles ⩽100°. Only events belonging to the Pass 6 "Diffuse" class are retained, as those events have the highest probability of being photons (Atwood et al. 2009). For the event probability weighting described below we use the P6_V11 Instrument Response Functions (IRFs). Using the gtmktime tool, times when the rocking angle of the satellite exceeded 52° are excluded. We also require that DATA_QUAL and LAT_CONFIG are unity and that Earth's limb does not impinge upon the region of interest (ROI).

As demonstrated in previous work (Bickel et al. 2008; Kerr 2011a), the sensitivity of gamma-ray pulsation searches can be improved by weighting the jth event with a probability w_j ∈ [0, 1] that it originated from a putative pulsar. The photon weights w_j are calculated by using a full spectral model of the region around the gamma-ray source, and by exploiting the IRFs to provide background rejection which is superior to simple angular and energy data-selection cuts. For the first time in this paper such a photon-weighting scheme has been applied in a blind search; further details are given in Section 3.

To calculate the weights w_j for each event from a given source, we perform likelihood spectral analyses using the pointlike tool (see Kerr 2011b, for a description). For each selected source, we construct a spectral model for the region by including all sources of the 2FGL catalog found within 8° of the selected source, using the spectral forms given in the catalog. The spectra of the selected sources are modeled as exponentially cutoff power laws, typical of known gamma-ray pulsars, of the form N₀(E/GeV)^−Γexp (− E/E_c), where N₀ is a normalization factor, Γ is the photon index, and E_c is the cutoff energy. The Galactic and extragalactic diffuse emission and residual instrument background also enter the calculation of the weights. The Galactic diffuse emission is modeled using the gll_iem_v02_P6_V11_DIFFUSE map cube, while the extragalactic diffuse and residual instrument backgrounds are modeled using the isotropic_iem_v02_P6_V11_DIFFUSE template (a detailed description of these background models can be found in Section 3 of Abdo et al. 2010a). These models are available for download at the Fermi Science Support Center.²⁵ The tool gtsrcprob is then used to calculate the event weights w_j based on the best-fit spectral models obtained from the maximum likelihood analyses.

In a year, the LAT detects of order 10³ photons from a typical gamma-ray pulsar; in the same year, a typical pulsar rotates at least 10⁸ times around its axis. The blind-search problem is to find a rotational-phase model $\Phi (t) = 2 \pi (f t + {\dot{f}} t^2 / 2)$ and a sky position that match the solar system barycenter (SSB) arrival times t of the different photons, where Φ denotes the rotation angle of the star about its axis, in radians, measured from its starting position at t = 0, and observed at the SSB. The signal hypothesis is that the photons' arrival times are "clustered" near specific "orientations" of the star (i.e., Φ(t) mod 2π deviates from uniformity on the interval [0, 2π]). The null hypothesis is that the arrival times of the photons are a random Poisson process. In this paper, we do not explicitly indicate the dependence of Φ on f, $\dot{f}$, and sky position, but this dependence is important and implicit in many formulae below.

To find a matching phase model, a grid of "templates" in the four-dimensional parameter space of sky position and $(f, \dot{f})$ is constructed. Note that the 2FGL catalog sky positions of the targeted unassociated sources based on the spatial distribution of events are typically not precise enough for pulsar searches. A search grid of sky points around this catalog position is needed to reduce signal loss arising from imperfect correction of the Doppler shifts caused by Earth's orbital motion around the SSB. The need for sky gridding is particularly acute for MSP spin frequencies. Therefore, in contrast to previously published blind searches,²⁶ we grid a circular sky region centered on the 2FGL catalog source location using a radius which is 20% larger than the semi-major axis of the 95% confidence elliptical error region (given by the "Conf_95_SemiMajor" parameter).

Unfortunately the number of templates (grid points) required to discretely cover the entire four-dimensional search parameter space increases as a high power of the coherent integration time. Hence a fully coherent approach for several years of data is computationally impossible. Therefore, we employ a search strategy which is designed to achieve maximum overall sensitivity at fixed computing cost.²⁷

To efficiently scan through years of Fermi-LAT data for previously unknown gamma-ray pulsars, we use a so-called hierarchical search approach. This is analogous to hierarchical methods used in searches for gravitational-wave pulsars (Schutz & Papa 1999; Papa et al. 2000; Brady & Creighton 2000; Abbott et al. 2009a, 2009b; Cutler et al. 2005). In a first semi-coherent stage, we here adopt the optimal metric-based gridding methods described in Pletsch & Allen (2009) along with the sliding coherence window technique (Pletsch 2011). In a second stage, significant semi-coherent candidates are automatically followed up in a fully coherent analysis. Finally, a third stage further refines coherent pulsar candidates by using higher harmonics. Full details of the complete search scheme will be presented in forthcoming work (H. J. Pletsch & L. Guillemot 2011, in preparation).

Here we first describe the principle of the method, to firmly establish the analogy with the existing gravitational-wave literature. Then we describe what is done in practice, which is mathematically equivalent (up to justifiable approximations) but computationally more efficient.

In the first stage, a semi-coherent detection statistic S is computed for each template. We refer to S as "semi-coherent" because it is effectively the incoherent sum over several years, of terms which are coherent over several days. The coherent terms are the power in Fourier bins, calculated by treating each photon arrival as a delta function in time.

Denoting the arrival time of the jth event (photon) at the SSB by t_j, the coherent power P_τ in a (Gaussian) window centered at time τ is defined by

$$\begin{equation} P_\tau = \biggl | \sum _{j=1}^N w_j \; e^{-i \Phi (t_j)} e^{-2 \pi (t_j - \tau)^2/T^2}\biggr |^2. \end{equation} \tag{ 1 }$$

The sum is taken over all photons (here, N = 8000) in the data set; the effective window duration is $\int e^{-2 \pi \tau ^2/T^2} d\tau = T/\sqrt{2}$. As described in Section 2, the weights w_j estimate the probability that the photon comes from the selected source.

To form the semi-coherent detection statistic S, the values of P_τ are summed ("incoherently combined")

$$\begin{equation} S = \frac{2}{T} \int d\tau \; P_\tau - \sum _{j=1}^N w_j^2. \end{equation} \tag{ 2 }$$

Note that in this definition we have subtracted a constant (phase model independent) term. Because it contains a Gaussian window, the integrand in Equation (2) falls off exponentially at early and late time (large values of |τ|). Thus, the limits of integration can be taken as the entire real line τ ∈ (− ∞, ∞); to good approximation this gives the same value as integrating only over the total observation interval (about 975 days in this search).

The semi-coherent detection statistic S is an incoherent sum of powers, which discards the phase information over time periods longer than of order T. This uniform overlap maximizes the search sensitivity for fixed T and for fixed computational resources (Pletsch 2011). For computational efficiency, in this search we choose the N = 8000 photons with the highest probabilities (largest values of w_j).

To understand how to compute S efficiently, one can explicitly evaluate Equation (2). Completing the squares in the product of the Gaussians and carrying out the integration over τ, one obtains

$$\begin{equation} S = \sum _{j=1}^N \sum _{k=1}^N w_j w_k {e}^{-i [\Phi (t_j) - \Phi (t_k)]} {e}^{-\pi (t_j - t_k)^2/T^2} - \sum _{j=1}^N w_j^2. \end{equation} \tag{ 3 }$$

Here the effective duration of the Gaussian window is $\int e^{-\pi \tau ^2/T^2} d\tau = T$. In practice, to compute the semi-coherent power efficiently, we replace the Gaussian window in Equation (3) with a rectangular window of the same duration T, as given below in Equation (5). In this search, the width of the rectangular window is T = 2¹⁹ s (≈6 days).

The template grid in parameter space is the Cartesian product of a rectangular two-dimensional grid in $(f, \dot{f})$ and a sky grid which has constant density when orthogonally projected onto the ecliptic plane. The problem of constructing efficient search grids has been intensively studied in the context of gravitational wave searches (see, e.g., Brady et al. 1998; Brady & Creighton 2000; Prix 2007; Pletsch & Allen 2009; Pletsch 2010) and we employ these concepts here. The values of frequency are equally spaced, separated by the fast Fourier transform (FFT) frequency-bin width Δf = 1/T. In practice, to reduce the fractional loss in S for frequencies not coinciding with Fourier frequencies, we use a computationally efficient interpolation, referred to as "interbinning" (Ransom et al. 2002). The spacing in the other three dimensions is determined by a metric which measures the fractional loss in the expected value of S that arises if the signal is not located exactly at a grid point (Balasubramanian et al. 1996; Owen 1996; Prix 2007). In spin-down $\dot{f}$, we use a uniform grid spacing $\Delta {\dot{f}} = {\sqrt{720 m}}/{\pi \gamma T^2}$. In this search m = 0.3 is the maximum tolerable fractional loss in S, and from Pletsch & Allen (2009),

$$\begin{equation} \gamma ^2 = 1 + \frac{60}{N} \sum _{j=1}^N \frac{(t_j-\bar{t})^2}{T^2}, \end{equation} \tag{ 4 }$$

where $\bar{t} = \sum _j t_j/N$ is the mean photon arrival time. The description of this grid can be found in Pletsch & Allen (2009) with a detailed derivation in Pletsch (2010). The grid in the sky is determined by the same metric, permitting a maximum fractional loss m in the value of S. The spacing of the sky grid is determined by the Doppler shift arising from Earth's (more precisely, the Fermi satellite's) motion around the Sun. At the north Ecliptic pole the angular spacing is $\Delta \theta = {\sqrt{2m}}\,c/{\pi f D}$, where c is the speed of light and D is a baseline distance (defined below). When the entire sky grid is projected into the plane of the ecliptic, the grid points are uniformly spaced on the plane (Astone et al. 2002; Abbott et al. 2009a, 2009b). This angular spacing is similar to the angular spacing in the diffraction pattern of a two-slit system, where the wavelength is c/f and the separation of the two slits is the straight-line distance D between two points on Earth's orbit about the Sun. If the coherent integration time T is less than half a year, then D = (998 s) c sin (πT/1 yr). If the coherent integration time T is greater than half a year, then D = (998 s) c is the diameter of Earth's orbit about the Sun.

To compute S efficiently, a time series is constructed and subsequently Fourier-transformed into the frequency domain. The time series contains T Δf_BW bins, where Δf_BW is the total frequency bandwidth being searched at a time using complex heterodyning²⁸ at the center of Δf_BW (see, e.g., Patel et al. 2010). The time series is initialized to zero, then the values of $w_j\, w_k\, e^{- i \pi \dot{f} \; (t^2_j-t^2_k) }$ are added into the bins determined by the time differences Δt_jk = t_j − t_k, for all pairs of photons j, k for which 0 < |Δt_jk| ⩽ T; the bin index is obtained by rounding the absolute value of the product Δf_BW Δt_jk to its nearest integer value.

Then the array is Fourier-transformed into the frequency domain (exploiting the FFT) to obtain S over the entire f grid. Up to an overall window-dependent normalization, one can write for given values of f, $\dot{f}$, and sky position,

$$\begin{equation} S = \sum _{j,k=1}^N Q(|\Delta t_{jk}|/T) \; w_j\, w_k\, e^{-2\pi i f\Delta t_{jk}- i \pi \dot{f} \; (t^2_j-t^2_k) }, \end{equation} \tag{ 5 }$$

where the rectangular function Q(x) is unity if 0 < x ⩽ 1 and vanishes otherwise.²⁹

Although other aspects are different, the use of an FFT applied to time differences is very similar to techniques previously used in blind searches of Fermi-LAT data (Atwood et al. 2006, Equation (3) therein has a typo which is corrected in Equation (2) of Ziegler et al. 2008). The method used in Atwood et al. (2006) was the first application of this classic method (e.g., Blackman & Tukey 1958) to gamma-ray astronomy (in estimating the power spectrum an approximate autocovariance function is calculated using a maximum lag, and then Fourier-transformed).

In contrast to previous searches, our method uses an optimal gridding of the parameter space for both the semi-coherent and coherent stages, as well as an automated follow-up, and incorporates the spin-down corrections in a way that permits heterodyning and highly efficient code.

The search was done on the 1680-node Atlas Computing Cluster (Aulbert & Fehrmann 2008) built around four-core processors with 8 GB of random-access memory; for these we used a heterodyning bandwidth Δf_BW = 256 Hz. Breaking the full frequency range of the search into frequency bands allows the computation to fit into memory and also allows the use of different sky grids in each band. This further reduced the computational cost, since the number of required sky grid points increases with the square of frequency.

After computing S on the four-dimensional grid in parameter space, points with statistically significant values of S are candidates for possible pulsar signals, and are followed up in a second stage. This is done by "refining the grid" and increasing the coherent integration time. This is a hierarchical scheme which is analogous to "zooming": successively swapping microscope objectives for ones of higher magnification, then recentering the interesting point on the slide (see, e.g., Cutler et al. 2005; Krishnan et al. 2004). In our case, this is done by constructing the fully coherent detection statistic P over the entire data set (or equivalently taking the Gaussian window-size T → ∞ in Equation (1)) obtaining

$$\begin{equation} P = \frac{1}{{\kappa }^2} \biggl | \sum _{j=1}^{N} w_j \; e^{-i \Phi (t_j)} \biggr |^2, \end{equation} \tag{ 6 }$$

where for convenience we have normalized by the positive constant κ given by

$$\begin{equation} {\kappa }^2 = \frac{1}{2} \sum _{j=1}^{N} w_j^2. \end{equation} \tag{ 7 }$$

The computing cost to coherently follow up a single candidate is negligible in comparison to the cost of the previous semi-coherent search.

In selecting statistically significant semi-coherent candidates which are automatically followed up in the second stage using a fully coherent analysis, we do not use a fixed threshold to define "statistical significance." In the semi-coherent stage, the search code keeps an internal list of the strongest signal candidates. Each member of this list is coherently followed up and corresponds to the largest value of S detected in eight adjacent spin-down values for the entire heterodyning frequency bandwidth and a single sky point.

The refined grid of the fully coherent follow-up covers a region of parameter space of size $(4 \Delta f) \times (4 \Delta {\dot{f}}) \times (4 \Delta \theta \times 4 \Delta \theta)$ when projected into the ecliptic plane. In other words, it covers a region whose volume is 256 times larger than the volume of a fundamental cell in the original grid: its extent in each dimension of parameter space is four grid intervals. The refined grid has a spacing given by the previous formulae for Δf, $\Delta {\dot{f}}$, and Δθ, except that the coherence time T is set equal to the length of the entire data set, and γ = 1. Since in this case only a small parameter-space region around the candidate is explored, it is computationally efficient to compute P directly in the time domain (FFTs are not used), exploiting the sparsity of the photon data.

If the value of P, which measures the fully coherent power (in a single harmonic), is statistically significant, then in a third stage further refinement is carried out using higher harmonics (Fourier components). We adopt the so-called H-test, which has been widely used in X-ray and gamma-ray pulsar detection (de Jager et al. 1989; de Jager & Büsching 2010). This test measures the statistical significance of the energy in the first 20 (non-DC) Fourier components of the pulse profile as a function of phase. As in Equation (1), the H-test can also be modified to include the photon probability weights w_j (see Kerr 2011a). Note that Equation (5) in Kerr (2011a) contains an error; corrected formulae used in this work are given below.

The weighted H-test statistic is defined as follows. For each photon arrival time t_j, the pulse phase x_j (between zero and one) is calculated as x_j = [Φ(t_j) mod 2π]/2π. The pulse profile (for 0 ⩽ x < 1) is a sum over the photons

$$\begin{equation} p(x) = \sum _{j=1}^N w_j \, \delta (x-x_j), \end{equation} \tag{ 8 }$$

where x_j is the pulse phase of the jth photon and δ(x) is a one-dimensional Dirac delta function. It can be expressed as a Fourier series

$$\begin{equation} p(x) = \kappa \sum _{\ell = -\infty }^\infty \alpha _\ell \, e^{2\pi i\ell x}, \end{equation} \tag{ 9 }$$

which implies that the (complex) Fourier coefficients α_ℓ are given by

$$\begin{equation} \alpha _\ell = \frac{1}{\kappa } \sum _{j=1}^N w_j \, e^{-2\pi i \ell x_j}, \end{equation} \tag{ 10 }$$

where the definition of κ is identical to Equation (7). The normalization of the Fourier coefficients has been chosen so that if the photon arrival times are uniformly distributed, independent random variables, then in the limit of the large numbers of photons, for ℓ > 0, ℜ(α_ℓ) and ℑ(α_ℓ) are independent Gaussian random variables with zero mean and unit variance. Finally, the weighted H-test statistic is defined as

$$\begin{equation} H = \max _{1 \le L \le 20} \biggl [\sum _{\ell =1}^L \bigl |\alpha _{\ell } \bigr |^2 - 4(L -1)\biggr ]. \end{equation} \tag{ 11 }$$

The quantity which is subtracted, 4(L − 1), is motivated by an empirical numerical study (de Jager et al. 1989), providing the best omnibus test for unknown pulse profiles.

Maximizing H over sky position, frequency f, and spin-down rate $\dot{f}$ selects the "narrowest and sharpest" overall pulse profile. In contrast, maximizing the statistic P = |α₁|² of Equation (6) favors putting more power into the lowest harmonics.

Using H as the test statistic, a further stage of parameter-space refinement is done in the same way as before: in the frequency, spin-down, and sky parameters, we cover regions that include four grid steps (in each dimension) of the previous grid. The chosen refinement factor (ratio of the number of grid points after and before refinement, in each dimension) is about an order of magnitude. At each grid point, the weighted H-statistic of Equation (11) is found. The parameter-space point with the largest statistic is selected as our best estimate of the pulsars' parameters which are then further refined through the timing-analysis procedure (Ray et al. 2011) described in Section 5.2.

The hierarchical search pipeline has been validated by successfully recovering previously known gamma-ray pulsars, including some of the brightest gamma-ray MSPs (Abdo et al. 2010c; H. J. Pletsch & L. Guillemot 2011, in preparation). In the next section, the complete search scheme is illustrated with a detailed example.

Before giving the results for all of the pulsars that have been discovered with this new search method, we first go through a single example in detail. This illustration uses PSR J1620–4927, the first new pulsar found in this work.

Figure 1(a) shows the first stage of the analysis. For each f and $\dot{f}$ value, the largest value of S found in the sky grid is displayed as an intensity. The point of highest intensity corresponds to PSR J1620–4927's f and $\dot{f}$ parameters.

Zoom In Zoom Out Reset image size

The second analysis stage is the automated follow-up of all candidates shown in Figure 1(a). As explained in Section 3, this is accomplished by carrying out a fully coherent search over a small region of parameter space around each candidate which has been identified as statistically significant in the previous semi-coherent stage. If the candidate found in the semi-coherent step was simply a statistical outlier, then the fully coherent statistic P will not be significant. Figure 1(b) presents the results of the fully coherent statistic P for PSR J1620–4927. Again the point of highest intensity compared to the background, indicating the presence of a coherent signal in the data set, is due to the new pulsar.

The importance of the photon probability weights w_j can be illustrated by repeating the analysis using the same 8000 photons with the weights set to unity. (Note that the weights were used in selecting the 8000 photons, so this is not a complete comparison.) The result is that Figure 1(a) still shows a statistically significant outlier in the semi-coherent search step, which is followed up automatically and gives a statistically significant outlier in the fully coherent output of Figure 1(b). In both steps the signal and its statistical significance are reduced, but the pulsar is still found. However for some of the other new pulsars, this is not the case: if the weights are set to unity, then the pulsar is not detected. Overall the weights play a larger role for sources in crowded regions of the sky, for example, near the Galactic plane.

At the third stage, further refinement is carried out by maximizing the weighted H-test statistic over a small region of parameter space around the most significant candidate from the previous step. As described in Section 3, the parameter-space grid used at this stage is yet another order-of-magnitude finer than the one of the previous stage. Figure 2 shows the corresponding results for PSR J1620–4927.

Zoom In Zoom Out Reset image size

The parameter-space point with the largest H-test statistic is selected as our best estimate of the pulsar parameters at this stage. The uncertainties in these estimated parameters can be obtained by using the fact that a 1σ deviation for a Gaussian distribution has a value ≈0.6 of the maximum. In the current example, as shown in Figure 2, the value of the maximum weighted H-test statistic is 566.

Figure 2(c) illustrates the way the weighted H-test statistic accumulates over the total observation time. The linearity of this plot indicates that the properties of the source, and of the data, were time stationary; it shows that the pulsar's emission is not changing with time. In this way it also gives additional confidence that the pulsar is real and not a statistical noise outlier. (We have not provided analogous plots for the other eight pulsars reported in this paper, but they are similar.)

As shown in de Jager & Büsching (2010) and Kerr (2011a), the weighted H-test has a false alarm probability approximately described by P_FA ≈ e^−0.4 H. For this example, the associated false alarm probability is log₁₀P_FA ≈ −98. We have not tried to rigorously estimate the trial factor, but since the total number of floating-point computations that could be performed by our computing systems over a period of some months is less than 10²¹, we conclude that PSR J1620–4927 is a real gamma-ray pulsar and not a statistical outlier.

Table 1 shows the names and sky positions of the nine discovered pulsars, and also lists known source associations. The inferred rotational parameters of the new pulsars are presented in Table 2. In addition, further derived parameters are given, including the spin-down power $\dot{E} = -4\pi ^2 I f \dot{f}$, where I = 10⁴⁵ g cm² is the assumed neutron star moment of inertia, and the magnetic field strength at the neutron star surface $B_{{\rm S}} \approx 3.2\times 10^{19} (-\dot{f}/f^3 \,\hbox{s}^{-1})^{1/2}$ G and at the light cylinder $B_{{\rm LC}}\approx 2.94\times 10^{8} (-\dot{f} f^3\,\hbox{s}^5)^{1/2}$ G, respectively.

Table 2. Measured and Derived Parameters of the Discovered Gamma-Ray Pulsars

Pulsar Name	f	$\dot{f}$	Weighted	τ	$\dot{E}$	BS	BLC
	(Hz)	(−10−13 Hz s−1)	H-test	(kyr)	(1034 erg s−1)	(1012 G)	(kG)
J0106+4855	12.02540173638(8)	0.61881(7)	843.1	3081.1	2.9	0.2	3.0
J0622+3749	3.00112633651(5)	2.28985(4)	288.8	207.8	2.7	2.9	0.7
J1620–4927	5.81616320951(5)	3.54782(4)	566.4	259.9	8.1	1.4	2.4
J1746–3239	5.01149235750(3)	1.64778(3)	249.8	482.2	3.3	1.2	1.3
J1803–2149	9.4044983174(2)	17.25894(6)	451.9	86.4	64.1	1.5	11.0
J2028+3332	5.65907208453(2)	1.55563(2)	1108.3	576.8	3.5	0.9	1.5
J2030+4415	4.4039248637(5)	1.2576(2)	584.8	555.2	2.2	1.2	1.0
J2111+4606	6.3359340865(4)	57.4218(3)	554.3	17.5	143.6	4.8	11.1
J2139+4716	3.5354509962(2)	0.2232(2)	351.1	2511.5	0.3	0.7	0.3

Notes. The reference epoch for all measured rotational parameters is MJD 55225 and the time range for all timing models is MJD 54682–55719. The derived quantities in Columns 5–8 are based on the f and $\dot{f}$ values obtained from the timing model and are rounded to the nearest significant digit. To model the timing noise present in PSR J1803–2149, a second frequency derivative is necessary: $\ddot{f}=$ 7.3(8) × 10⁻²⁴ Hz s⁻². To model the timing noise present in PSR J2030+4415, higher frequency derivatives up to third order are necessary: $\ddot{f}=-$1.5(3) × 10⁻²³ Hz s⁻² and $\dddot{f}=-$6(2) × 10⁻³¹ Hz s⁻³. To model the timing noise present in PSR J2111+4606, higher frequency derivatives up to fourth order are necessary: $\ddot{f}=$ 2.30(5) × 10⁻²² Hz s⁻², $\dddot{f}=-$7.9(2) × 10⁻³⁰ Hz s⁻³, and $\ddddot{f}=$ 3.2(4) × 10⁻³⁷ Hz s⁻⁴. The numbers in parentheses are the statistical 1σ errors in the last digits.

Download table as:  ASCIITypeset image

Figure 3 plots the newly discovered pulsars on an f–$\dot{f}$ diagram, where they can be compared with the known pulsar population. One can see that the new pulsars lie in a similar region (and hence belong at large to the same population) as the previously found blind-search Fermi-LAT gamma-ray pulsars.

Zoom In Zoom Out Reset image size

Two of the nine new systems, PSRs J1803–2149 and J2111+4606, are young energetic pulsars ($\dot{E} \gtrsim 6\times 10^{35}$ erg s⁻¹ and τ ≲ 100 kyr), located near the Galactic plane. Among the seven remaining less-energetic and older pulsars, five are also located near the Galactic plane (|b| < 5°) and two are found at higher Galactic latitudes (|b| > 10°). One of these, PSR J0106+4855, has the largest characteristic age τ = 3 Myr and the lowest surface magnetic field strength (B_S ≈ 2 × 10¹¹ G) of all LAT blind-search pulsars found to date. Also standing out, PSR J2139+4716 has the smallest spin-down power ($\dot{E} = 3\times 10^{33}$ erg s⁻¹) among all non-recycled gamma-ray pulsars ever detected (see De Luca et al. 2011).

5.1. Source Associations

5.2. Timing Analysis

For each new gamma-ray pulsar, the definitive parameters as given in Tables 1 and 2 are determined via precise measurements of pulse times of arrival (TOAs) and by fitting the parameters of a timing model to these measurements using methods described in Ray et al. (2011).

To obtain a set of TOAs, photons are extracted for each source using a radius and minimum energy cut to optimize the signal-to-noise ratio for each pulsar. Then these data are subdivided into segments of about equal length. The set of best-guess pulsar parameters from the search (maximizing H) is used as a preliminary ephemeris to fold the photon arrival times, producing a set of pulse profiles. The TOAs are then measured by cross-correlating each pulse profile with a multi-Gaussian or kernel-density template derived from fitting the entire data set. This is done using the unbinned maximum likelihood method described in Ray et al. (2011). Then Tempo2 (Hobbs et al. 2006) is used to fit the TOAs to a timing model including sky position, frequency, and frequency derivative.

As often is the case, the youngest object found in this work (PSR J2111+4606, τ = 17 kyr) exhibits particularly large timing noise. In turn, this requires including higher-order frequency-derivative terms to make a good timing model fit; in this case up to the fourth derivative of the frequency as shown in Table 2. For the same reason, the timing models in Table 2 for PSRs J2030+4415 and J1803–2149 include terms up to the third and the second frequency derivative, respectively.

Based on these timing solutions, Figures 4–12 show the resulting phase–time diagrams and pulse profiles for each of the newly discovered pulsars. These plots are obtained from calculating the phase for each of the 8000 photons selected for the blind search and using the parameters listed in Tables 1 and 2, and weighting each event by its probability of originating from the pulsar w_j. The integrated pulse profiles (weighted pulse phase histograms) are constructed with a resolution of 32 bins in phase per rotation. The 1σ error bars in the integrated pulse profiles are statistical only and are given by (∑_jw²_j)^1/2, where j runs over all events falling into the same phase bin. (Note that the formula for the fractional error has the opposite sign of the exponential: −1/2.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As seen in Figures 4–12, eight of the nine pulsars have pulse profiles with two peaks. For pulsars where the two peaks are separated by nearly one-half of a rotation, it is possible to detect or discover the pulsar at the second harmonic (i.e., at twice the actual spin frequency of the pulsar). For pulsar candidates which were discovered with (apparently) single-peaked profiles, we tested for the true fundamental spin frequency by folding at the subharmonic of the discovery frequency. If the subharmonic is the correct frequency, the two resulting peaks may satisfy one or more of the following conditions: offsets that are measurably different than 0.5 in phase, significant differences in the integrated weighted counts under each peak, or significant differences in the shape of the peaks in different energy bands. With modest signal-to-noise ratios, such determinations are not always conclusive. For PSR J0106+4855, our identification of the true period was subsequently confirmed by the detection of radio pulsations with a single-peaked profile. For the one apparently single-peaked pulsar in our sample, PSR J2139+4716, none of the above tests yield strong evidence for the profile being double peaked at half the frequency. Additional data will be required to strengthen this conclusion.

In order to further characterize each pulse profile, we fit the pulsars' weighted gamma-ray peaks to Lorentzian lines. The derived pulse shape parameters are listed in Table 3. Note that PSRs J0622+3749 and J1746–3239 show indications of sub-structures in their main gamma-ray peak. For these pulsars a single Lorentzian line function is used for fitting the main component. Apart from PSR J0106+4855, which is detected at radio wavelengths (see Figure 13), there is no particular reference for absolute phase of these pulsars. For the other eight pulsars we have arbitrarily assigned the absolute phase reference such that the first peak occurs at a value of 0.1 in phase. The gamma-ray pulse profiles shown in Figures 4–12 along with the pulse-profile parameters listed in Table 3 are very similar to those of the previously discovered gamma-ray pulsars in blind searches (Abdo et al. 2009a; Saz Parkinson et al. 2010; Abdo et al. 2010c), further supporting the theory that the gamma-ray emission consists of fan beams produced in the outer magnetosphere.

Zoom In Zoom Out Reset image size

Table 3. Pulse-Profile Parameters of the Discovered Gamma-Ray Pulsars

Pulsar Name	Peak Multiplicity	FWHM1	FWHM2	Δ
J0106+4855	2	0.02 ± 0.01	0.05 ± 0.01	0.50 ± 0.01
J0622+3749	2	0.21 ± 0.04	0.03 ± 0.03	0.47 ± 0.03
J1620–4927	2	0.13 ± 0.04	0.15 ± 0.05	0.21 ± 0.03
J1746–3239	2	0.32 ± 0.07	0.11 ± 0.11	0.63 ± 0.06
J1803–2149	2	0.04 ± 0.02	0.05 ± 0.03	0.40 ± 0.02
J2028+3332	2	0.10 ± 0.02	0.04 ± 0.01	0.38 ± 0.02
J2030+4415	2	0.09 ± 0.03	0.04 ± 0.01	0.49 ± 0.02
J2111+4606	2	0.08 ± 0.02	0.07 ± 0.02	0.31 ± 0.02
J2139+4716	1	0.13 ± 0.03	...	...

Notes. For each of the nine pulsars, we give the parameters describing the shape of the pulse profile, including the peak multiplicity, the full widths at half maxima (FWHM) of the peaks, and the separation Δ between the gamma-ray peaks for pulsars with more than one peak.

Download table as:  ASCIITypeset image

5.3. Spectral Parameters

The spectral parameters for the new pulsars are obtained by fitting each phase-averaged spectrum with an exponentially cutoff power law with a photon index Γ and a cutoff energy E_c. The results for each pulsar are listed in Table 4. In addition to Γ and E_c which are explicit parameters of the fit, Table 4 also gives the important derived physical quantities of photon flux F₁₀₀ (in units of photons cm⁻² s⁻¹) and the energy flux G₁₀₀ (in units of erg cm⁻² s⁻¹) for events with energies between 100 MeV and 100 GeV.

Table 4. Spectral Parameters of the Discovered Gamma-Ray Pulsars

Pulsar Name	Γ	Ec	F100a	G100b	Lps	dps
		(GeV)	(10−8 photons cm−2 s−1)	(10−11 erg cm−2 s−1)	(1033 erg s−1)	(kpc)
J0106+4855	1.47 ± 0.23 ± 0.12	3.31 ± 0.92 ± 0.08	2.56 ± 0.77 ± 0.32	2.40 ± 0.31 ± 0.04	5.5	1.4c
J0622+3749	0.59 ± 0.34 ± 0.09	0.60 ± 0.13 ± 0.04	2.21 ± 0.35 ± 0.08	1.69 ± 0.15 ± 0.05	5.3	1.6
J1620–4927	1.01 ± 0.18 ± 0.05	2.44 ± 0.42 ± 0.28	9.61 ± 1.68 ± 0.94	13.5 ± 1.0 ± 1.7	9.1	0.7
J1746–3239	1.33 ± 0.08 ± 0.35	1.65 ± 0.12 ± 0.52	9.97 ± 0.94 ± 1.85	7.86 ± 0.41 ± 0.77	5.8	0.8
J1803–2149	1.96 ± 0.11 ± 0.20	5.73 ± 1.72 ± 2.07	20.7 ± 3.1 ± 0.5	13.1 ± 1.1 ± 2.1	25.6	1.3
J2028+3332	0.86 ± 0.21 ± 0.07	1.53 ± 0.24 ± 0.08	5.12 ± 0.87 ± 0.44	6.09 ± 0.41 ± 0.13	6.0	0.9
J2030+4415	1.89 ± 0.14 ± 0.22	2.16 ± 0.65 ± 0.67	13.3 ± 1.4 ± 0.2	7.06 ± 0.48 ± 0.66	4.8	0.7
J2111+4606	1.63 ± 0.14 ± 0.05	5.43 ± 1.80 ± 1.56	4.39 ± 0.69 ± 0.02	4.13 ± 0.34 ± 0.30	38.4	2.7
J2139+4716	0.80 ± 0.27 ± 0.02	1.02 ± 0.21 ± 0.07	2.65 ± 0.44 ± 0.19	2.51 ± 0.21 ± 0.01	1.8	0.8

Notes. This table describes the spectral properties of each of the nine pulsars, modeling each spectrum as an exponentially cutoff power law with photon indices Γ and cutoff energies E_c. The spectral parameters listed here for each pulsar are obtained from maximum likelihood fits. The first quoted uncertainties are statistical, while the second are systematic and correspond to the differences in the best-fit parameters observed when doing the spectral analyses with the P6_V3 IRFs and associated diffuse emission models (namely, the gll_iem_v02 map cube and isotropic_iem_v02 template). For each object, the pseudo-gamma-ray luminosity L_ps and the pseudo-distance d_ps are inferred from the apparent spin-down power $\dot{E}$ and the energy flux G₁₀₀ above 100 MeV. Note that these estimated gamma-ray luminosities and distances are subject to a number of caveats, detailed in Saz Parkinson et al. (2010), and could differ significantly from the actual values. ^aPhoton flux measured above 100 MeV. ^bEnergy flux measured above 100 MeV. ^cThe actual distance is 3.0 kpc, as inferred from the dispersion of the radio pulse measuring the free electron column density; see Section 6.1.

Download table as:  ASCIITypeset image

Analogous to the pulse-profile properties, Γ and E_c measured for the new pulsars are also similar to those observed for previously detected gamma-ray pulsars (Abdo et al. 2010c). This is not surprising, because (as described in Section 2) target sources for our search have been selected based on similarity of their spectral properties to known gamma-ray pulsars.

The distance (3 kpc) for one of the new pulsars (PSR J0106+4855) can be inferred based on the dispersion of the radio pulse measuring the free electron column density (see Section 6.1 for details). As the remaining eight pulsars are radio-quiet, this method cannot be used to estimate their distance. Furthermore, none of the pulsars are associated with a known SNR, preventing us from deriving distance estimates from such source associations.

However, it is still possible to obtain a crude estimate of the distance to the new pulsars, by exploiting the observed correlation between the gamma-ray luminosity L_γ and the spin-down power $\dot{E}$ for other gamma-ray pulsars with distance measures (see Abdo et al. 2010c). Based on this correlation "pseudo-gamma-ray luminosities" L_ps are derived as

$$\begin{equation} L_{\rm ps} \sim 3.2 \times 10^{33} ({\dot{E} / 10^{34}\, \rm {erg \ s}^{-1} })^{1/2}\, \rm {erg \ s}^{-1}, \end{equation} \tag{ 12 }$$

where the $\dot{E}$ values are obtained from Table 2. Assuming a geometrical correction factor $f_\Omega =1$ for the emission cone (Watters et al. 2009) for all gamma-ray pulsars, the relation $L_\gamma = 4 \pi f_\Omega \, G_{100}\, d^2$ is used to convert the energy flux and pseudo-gamma-ray luminosity into a "pseudo-distance" d_ps, following Equation (2) of Saz Parkinson et al. (2010):

$$\begin{equation} d_{\rm ps} \sim 1.6 \, \frac{({\dot{E} / 10^{34}\, \rm {erg \ s}^{-1} })^{1/4}}{(G_{100} / 10^{-11}\ \rm {erg} \ \rm {cm}^{-2} \ \rm {s}^{-1})^{1/2}} \, {\rm kpc}. \end{equation} \tag{ 13 }$$

For each of the new pulsars the resulting values for L_ps and d_ps are shown in Table 4. Note that these estimated gamma-ray luminosities and distances are subject to a number of caveats, detailed in Saz Parkinson et al. (2010), and could differ significantly from the actual values.

5.4. Why Were the New Pulsars not Found in Previous Blind Searches?

6.1. Pulsed Radio Emission

These discoveries represent a substantial increase in the number of gamma-ray pulsars detected in blind searches of Fermi-LAT data. Of the 26 previously discovered LAT blind-search pulsars, only 3 have been found to pulse in the radio band (Camilo et al. 2009; Abdo et al. 2010b; Saz Parkinson et al. 2010; P. M. Saz Parkinson et al. 2011, in preparation) and there are tight upper limits on the others (Ray et al. 2011; Keith et al. 2011). In order to exploit these new discoveries in population studies of the relative beaming fraction and geometry of the radio and gamma-ray emission, it is essential to determine if they are also visible as radio pulsars.

Because the source list was chosen from Fermi-LAT pulsar-like unassociated sources, all of these sources had been previously searched for radio pulsations by the Fermi Pulsar Search Consortium (Ransom et al. 2011; Hessels et al. 2011). We have re-analyzed all of these archival observations, with increased sensitivity because we can now do a single frequency trial folding the radio data with the gamma-ray ephemeris, and search only in a single parameter, the dispersion measure (DM). This greatly reduces the number of points searched in parameter space (the "trial factor") and implies that much smaller pulsed signal amplitudes are statistically significant. Where we saw an opportunity to go significantly deeper, we also made a number of new radio observations.

The telescopes and observing configurations used are described in Table 5, and the characteristics of the individual radio observations are shown in Table 6. We compute the sensitivities using the modified radiometer equation given in Equation (A1.22) of Lorimer & Kramer (2005):

$$\begin{equation} \mathcal {S}_\mathrm{min} = \beta \frac{ 5\, T_\mathrm{sys}}{G \sqrt{n_\mathrm{p} t_\mathrm{int} \Delta F}} \sqrt{\frac{W}{1/f-W}} , \end{equation} \tag{ 14 }$$

where β is the instrument-dependent factor due to digitization and other effects (when unknown, we assume β = 1.25), a value of 5 has been assumed for the threshold signal-to-noise ratio for a confident pulsar detection, T_sys = T_rec + T_sky, G is the telescope gain, n_p is the number of polarizations used (2 in all cases), t_int is the integration time, ΔF is the observation bandwidth, f is the pulsar spin frequency, and W is the pulse width (for uniformity, we assume W = 0.1/f). Note that conventionally the 3 K background temperature is included in the receiver temperature T_rec, which is measured on cold sky, and so T_sky represents the excess temperature from the Galactic synchrotron component, which we estimate by scaling the 408 MHz map of Haslam et al. (1982) to the observing frequency ν with a spectral index α = −2.6 (defined as $\mathcal {S}_\nu \propto \nu ^{\alpha }$).

Table 5. Definition of Radio Observing Codes

Obs Code	Telescope	Gain	Frequency	Bandwidth ΔF	βa	np	HWHM	Trec
		(K Jy−1)	(MHz)	(MHz)			(arcmin)	(K)
GBT-350	GBT	2.0	350	100	1.05	2	18.5	46
GBT-820	GBT	2.0	820	200	1.05	2	7.9	29
GBT-S	GBT	1.9	2000	700b	1.05	2	3.1	22
Eff-L1	Effelsberg	1.5	1400	250	1.05	2	9.1	22
Eff-L2	Effelsberg	1.5	1400	140	1.05	2	9.1	22
Jodrell	Lovell	0.9	1520	200	1.05	2	6.0	24
AO-327	Arecibo	11	327	25	1.12	2	6.3	116
AO-Lwide	Arecibo	10	1510	300	1.12	2	1.5	27
Parkes-BPSR	Parkes	0.735	1352	340	1.05	2	7.0	25

Notes. The sky locations of all nine pulsars have been searched for pulsating radio emissions. This table gives the radio telescope and back-end parameters used in those observations, which are described in Table 6. ^aInstrument-dependent sensitivity degradation factor. ^bThe instrument records 800 MHz of bandwidth, but to account for a notch filter for radio frequency interference and the lower sensitivity near the band edges, we use an effective bandwidth of 700 MHz for the sensitivity calculations.

Download table as:  ASCIITypeset image

Table 6. Radio Search Observations of the New Gamma-Ray Pulsars

Target	Obs Code	Date	tint	R.A.a	Decl.a	Offset	Tsky	$\mathcal {S}_\mathrm{min}$
			(minutes)	(J2000)	(J2000)	(arcmin)	(K)	(μJy)
J0106+4855	GBT-350	2009 Oct 25	32	01:06:37.7	48:54:11	2.7	49.3	136
	GBT-820	2010 Nov 17	45	01:06:35.5	48:55:30	1.8	5.4	30b
	GBT-820	2010 Dec 17	45	01:06:35.5	48:55:30	1.8	5.4	30b
	Eff-L2	2011 Jun 1	45	01:06:25.1	48:55:52	0.0	1.3	31
J0622+3749	GBT-350	2009 Oct 27	32	06:22:05.5	37:51:07	2.2	46.2	131
	Eff-L1	2010 May 14	32	06:22:14.7	37:51:49	2.8	1.3	30
	Eff-L1	2010 Feb 6	10	06:22:15.0	37:51:48	2.8	1.3	53
	Eff-L1	2010 Jul 10	52	06:22:13.0	37:50:36	1.5	1.3	22
	Eff-L1	2010 Jul 10	55	06:22:13.0	37:50:36	1.5	1.3	22
	GBT-820	2010 Dec 12	45	06:21:59.0	37:51:36	3.3	5.0	32
	GBT-820	2010 Dec 17	45	06:21:59.0	37:51:36	3.3	5.0	32
J1620−4927	Parkes-BPSR	2009 Aug 3	270	16:21:05.5	−49:30:32	4.9	16.9	42
	Parkes-BPSR	2010 Nov 18	144	16:20:43.5	−49:28:24	0.9	16.9	42
	Parkes-BPSR	2011 May 10	72	16:20:41.3	−49:27:36	0.0	16.9	58
J1746−3239	GBT-S	2009 Dec 23	60	17:46:47.9	−32:36:22	3.8	5.1	30
	GBT-820	2010 Nov 14	45	17:46:41.0	−32:36:18	4.6	51.6	85
	Parkes-BPSR	2011 May 10	72	17:46:54.9	−32:39:55	0.0	14.1	54
J1803−2149	Eff-L1	2010 Feb 13	25	18:03:12.0	−21:47:27	1.6	17.8	55
	Eff-L1	2010 May 22	32	18:03:11.7	−21:47:28	1.5	17.8	48
	GBT-S	2010 Sep 4	65	18:03:11.7	−21:47:28	1.5	7.1	13
J2028+3332	GBT-820	2009 Aug 13	60	20:27:48.0	33:32:24	6.6	10.3	46
	GBT-S	2010 Sep 20	30	20:28:18.0	33:33:23	1.4	1.0	15
	GBT-820	2010 Nov 22	45	20:28:19.0	33:32:53	0.8	10.3	33
	GBT-820	2010 Dec 17	45	20:28:19.0	33:32:53	0.8	10.3	33
	AO-Lwide	2011 May 21	45	20:28:19.9	33:32:06	0.0	2.1	4
	AO-327	2011 May 30	25	20:28:19.9	33:32:06	0.0	112.9	142
J2030+4415	Eff-L1	2010 Feb 7	10	20:30:55.0	44:11:52	3.8	4.0	62
	Eff-L1	2010 May 15	32	20:30:55.3	44:11:53	3.8	4.0	35
	Eff-L1	2010 Jul 10	60	20:30:59.2	44:15:33	1.4	4.0	23
	Eff-L1	2010 Jul 30	60	20:30:59.2	44:15:33	1.4	4.0	23
	GBT-820	2010 Nov 22	45	20:30:54.7	44:16:08	0.8	16.0	38
	GBT-820	2011 May 28	183	20:30:51.3	44:15:38	0.0	16.0	19
	Eff-L2	2011 Jun 1	45	20:30:51.5	44:15:37	0.0	4.0	35
J2111+4606	GBT-820	2009 Sep 19	60	21:11:22.8	46:05:53	0.3	16.0	33
	Jodrell	2011 Jun 22c	14 × 60	21:11:24.0	46:06:29	0.0	9.0	14
J2139+4716	Eff-L1	2010 Jul 10	60	21:39:52.3	47:13:43	2.6	2.0	22
	GBT-350	2009 Oct 25	32	21:39:53.2	47:15:22	1.0	74.9	171
	GBT-820	2010 Dec 11	45	21:39:53.5	47:13:30	2.7	8.2	34
	GBT-820	2010 Dec 18	45	21:39:53.5	47:13:30	2.7	8.2	34
	Eff-L1	2010 May 15	32	21:39:56.9	47:15:28	0.8	2.0	29

Notes. The sky locations of all nine pulsars have been searched for radio pulsations. Only for PSR J0106+4855 radio pulsations are detected. The minimum detectable flux density $\mathcal {S}_\mathrm{min}$ for each observation is computed at the observing frequency using Equation (14) and the parameters in Table 5, as described in the text. ^aTelescope pointing direction (not necessarily source position). ^bFor these two observations, radio pulsations were detected at a flux density of 20 μJy (see the text). ^cObserved 14 times for 1 hr each between this date and 2011 July 11.

Download table as:  ASCIITypeset image

We use a simple approximation of a telescope beam response to adjust the flux sensitivity in cases where the pointing direction was offset from the true direction to the pulsar. This factor is given by $q = e^{ -(\theta /\mathrm{HWHM})^2/1.5 }$, where θ is the offset from the beam center and HWHM is the beam half-width at half-maximum. A computed flux density limit of $\mathcal {S}_\mathrm{min}$ at the beam center is thus corrected via division by q for targets offset from the pointing direction.

For eight of the nine pulsars, we have established that they are indeed radio-quiet (or extremely radio-faint), as viewed from Earth.

For PSR J0106+4855 our observations and analysis have revealed very faint radio pulsations. In two of our archival 45 minute Green Bank Telescope (GBT) observations at 820 MHz, we detect the pulsar at a DM of 70.87 ± 0.2 pc cm⁻³ with a single-trial significance of 4σ–5σ in each observation. Although the detections are not individually very strong, we gain additional confidence in their veracity from the fact that the peak phases are consistent to 0.01 pulse periods when each observation is folded using the ephemeris determined from our gamma-ray timing. The summed radio profile shows a narrow peak with a duty cycle of ∼2%, as shown in Figure 13. We estimate a flux density of ∼20 μJy in both observations, using the standard radiometer equation and a measurement of the off-pulse noise level. These radio pulsations are detected at a flux density below the nominal detection limit of 30 μJy because the duty cycle is a factor of five smaller than the 10% duty cycle used in the sensitivity calculation for unknown pulse shapes. This corresponds to an equivalent flux density of 8 μJy at 1400 MHz using a typical pulsar spectral index α = −1.7. As seen in Table 6, we have made two other observations of this source, one at 350 MHz with the GBT and one at 1.4 GHz with Effelsberg. Neither of these observations detect the radio pulsations. Accounting for this narrow duty cycle of the pulse, the minimum detectable flux density for the Effelsberg observation at 1.4 GHz was about 15 μJy, so the non-detection in that observation is not surprising. This constrains the spectral index to be α ⩽ −0.5. On the other hand, a spectral index of −1.7 would imply a flux density at 350 MHz of 85 μJy, which is above the nominal sensitivity of our 350 MHz observation. This implies that either α > −1.3 or the sensitivity of that observation was affected by scatter or DM broadening or higher than expected sky background. Using the NE2001 model (Cordes & Lazio 2002), the measured DM corresponds to an estimated distance of 3.0 kpc, over a factor of two larger than the pseudo-distance. Given the radio detection, we can measure the so-called phase lag δ = 0.073 ± 0.003 between the gamma-ray and radio emission.

6.2. Continuum Radio Emission

6.3. X-Ray

This work reports on the discovery of nine gamma-ray pulsars through the application of a new blind-search method to about 975 days of Fermi-LAT data. The new pulsars were found by searching unassociated sources with typical pulsar-like properties selected from the Fermi-LAT Second Source Catalog (Abdo et al. 2011b).

The sensitivity of blind searches for previously unknown gamma-ray pulsars is limited by finite computational resources. Thus, efficient search strategies are necessary maximizing the overall search sensitivity at fixed computing cost. We have developed a novel hierarchical search method, adapted from an optimal semi-coherent method (Pletsch & Allen 2009) together with a sliding coherence window technique (Pletsch 2011) originally developed for detection of continuous gravitational-wave signals from rapidly spinning isolated neutron stars. The first stage of the method is semi-coherent, because coherent power computed using a window of six days is incoherently combined by sliding the window over the entire data set. In a subsequent stage, significant semi-coherent candidates are automatically followed up via a more sensitive fully coherent analysis. The method extends the pioneering methods first described in Atwood et al. (2006).

The new method is designed to find isolated pulsars up to kHz spin frequencies, by scanning a template grid in the four-dimensional parameter space of frequency, frequency derivative, and sky location. The optimal and adaptive sky gridding (not done in previously published blind searches) is necessary, particularly at the higher search frequencies, because the source-catalog sky positions are not precise enough to retain most of the signal-to-noise ratio in year-long data sets. A fundamental new element of the method is the exploitation of a parameter-space metric (well studied in the continuous gravitational-wave context; see, e.g., Brady & Creighton 2000; Prix 2007) to build an efficient template grid, as well as a metric approach to construct an optimal semi-coherent combination step (Pletsch & Allen 2009; Pletsch 2010). A further enhancement over previous searches is the sub-division of the total search frequency range into bands via complex heterodyning. This accommodates memory limitations, parallelizes the computational work, and permits the use of efficient sky grids adapted to the highest frequency searched in each band. A photon probability weighting scheme (Kerr 2011a) is also used for the first time in a published blind search, further improving the search sensitivity.

The nine newly discovered pulsars increase the previously known Fermi-LAT blind-search pulsar population (Abdo et al. 2009a; Saz Parkinson et al. 2010; P. M. Saz Parkinson et al. 2011, in preparation) by more than one-third and bring the total number to 35. The inferred parameters of the new pulsars suggest that they belong to the same general population as the previously found blind-search gamma-ray pulsars. Deep follow-up observations with radio telescopes have been conducted for all of the new pulsars, but significant radio pulsations have only been found for PSR J0106+4855. The null results for the other eight pulsars indicate that they belong to the growing population of radio-quiet gamma-ray pulsars, which can only be detected via their gamma-ray pulsations.

The computational work of the search has been done on the 6720-CPU-core Atlas Computing Cluster (Aulbert & Fehrmann 2008) at the Albert Einstein Institute in Hannover. Recently, in 2011 August, we have moved the computational burden of the search onto the volunteer distributed computing system Einstein@Home³³ (Abbott et al. 2009a, 2009b; Knispel et al. 2010). This will provide significantly more computing power and will allow a complete search of the parameter space up to kHz pulsar spin frequencies. We also hope that in the future an improved version of these methods can be used to carry out blind searches for gamma-ray pulsars in binary systems.

The combination of improved search techniques and much more powerful computational resources leaves us optimistic that we can find still more gamma-ray pulsars in the Fermi-LAT data. These advances should also greatly increase the chance of finding the first radio-quiet gamma-ray MSP with the Fermi-LAT. We hope that further discoveries, and further study of these systems, will eventually provide important advances in our understanding of pulsars, and of their emission mechanisms and geometry.

We thank Mallory Roberts, Ramesh Karuppusamy, Joris Verbiest, and Kejia Lee for help in retrieving and analyzing archival radio data, and for comments on and corrections to the manuscript, and to J. Eric Grove for similar assistance with archival X-ray data. We are grateful to Robert P. Johnson for his helpful comments on the manuscript, and to David Thompson for checking the EGRET 3EG catalogs, and for his management and organizational efforts on our behalf.

This work was partly supported by the Max Planck Gesellschaft and by U.S. National Science Foundation Grants 0555655, 0970074, and 1104902.

The National Radio Astronomy Observatory is 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 operation phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.