{'fraction_non_alphanumeric': 0.03935291770267366, 'fraction_numerical': 0.03496331767811646, 'mean_word_length': 4.4527575529333, 'pattern_counts': {'":': 1, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discovery of Nine Gamma-Ray Pulsars in Fermi-LAT Data Using a New Blind Search Method

We report the discovery of nine previously unknown gamma-ray pulsars in a blind search of data from the Fermi Large Area Telescope (LAT). The pulsars were found with a novel hierarchical search method originally developed for detecting continuous gravitational waves from rapidly rotating neutron stars. Designed to find isolated pulsars spinning at up to kHz frequencies, the new method is computationally efficient, and incorporates several advances, including a metric-based gridding of the search parameter space (frequency, frequency derivative and sky location) and the use of photon probability weights. The nine pulsars have spin frequencies between 3 and 12 Hz, and characteristic ages ranging from 17 kyr to 3 Myr. Two of them, PSRs J1803-2149 and J2111+4606, are young and energetic Galactic-plane pulsars (spin-down power above 6e35 erg/s and ages below 100 kyr). The seven remaining pulsars, PSRs J0106+4855, J0622+3749, J1620-4927, J1746-3239, J2028+3332, J2030+4415, J2139+4716, are older and less energetic; two of them are located at higher Galactic latitudes (|b|>10 deg). PSR J0106+4855 has the largest characteristic age (3 Myr) and the smallest surface magnetic field (2e11 G) of all LAT blind-search pulsars. PSR J2139+4716 has the lowest spin-down power (3e33 erg/s) among all non-recycled gamma-ray pulsars ever found. Despite extensive multi-frequency observations, only PSR J0106+4855 has detectable pulsations in the radio band. The other eight pulsars belong to the increasing population of radio-quiet gamma-ray pulsars.


INTRODUCTION
The Fermi Gamma-ray Space Telescope has been operating since launch in 2008 June. The Large Area Telescope (LAT) on board the Fermi satellite has an effective area of ∼ 0.8 m 2 (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 ∼ 0. • 8 ×(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 radiofrequency 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, 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. This paper describes the "blind" discovery of 9 additional pulsars using a new method. The blind search problem is computationally demanding because the relevant pulsar parameters (typically skyposition, frequency f , and spin-down rateḟ ) 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 "timedifferencing technique" as described in Atwood et al. (2006); 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 −13 Hz s −1 ≤ḟ ≤ 0, with a corresponding range of characteristic ages τ = − f /2ḟ . 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.

SOURCE SELECTION AND DATA PREPARATION
The Fermi-LAT Second Source Catalog (2FGL, Abdo et al. 2011b) lists 1873 sources, described by fits to ellipticallyshaped 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. Previouslyobserved 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). Gammaray pulsars tend to have large curvature significances ("Sig-nif_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) 25 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 tions (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 the 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 cut-off power laws, typical of known gamma-ray pulsars, of the form N 0 E/GeV −Γ exp −E/E c , where N 0 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 is 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 26 . 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.

THE NEW SEARCH METHOD
In a year, the LAT detects of order 10 3 photons from a typical gamma-ray pulsar; in the same year, a typical pulsar rotates at least 10 8 times around its axis. The blind-search problem is to find a rotational-phase model Φ(t) = 2π( f t +ḟ 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 ,ḟ 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 ,ḟ ) 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 26 http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html position is needed to reduce signal loss arising from imperfect correction of the Doppler shifts caused by the 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 27 , 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 28 .
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,b;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 (Pletsch & Guillemot 2011).
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 The sum is taken over all photons (here, N = 8000) in the data set; the effective window duration is e −2πτ 2 /T 2 dτ = T / √ 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 27 Ziegler et al. (2008) argues correctly that for two-week data-stretches, sky gridding is not essential. However for data stretches of length comparable to a year or longer, sky gridding is necessary to avoid significant loss of signal-to-noise ratio. 28 If this constraint is removed, then it is obvious that a more sensitive method exists. of P τ are summed ("incoherently combined") Note that in this definition we have subtracted a constant (phase model indepdendent) 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 Here the effective duration of the Gaussian window is e −πτ 2 /T 2 dτ = 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 19 s (≈ 6 days).
The template grid in parameter space is the Cartesian product of a rectangular two-dimensional grid in ( 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 FFT frequencybin 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ḟ , we use a uniform grid spacing ∆ḟ = √ 720m/πγT 2 . In this search m = 0.3 is the maximum tolerable fractional loss in S, and from Pletsch & Allen (2009), wheret = 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 maxi-mum fractional loss m in the value of S. The spacing of the sky grid is determined by the Doppler-shift arising from the Earth's (more precisely, the Fermi satellite's) motion around the Sun. At the north Ecliptic pole the angular spacing is ∆θ = √ 2m c/π 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,b). 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 the 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 the 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 29 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πḟ (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 ,ḟ and sky position, where the rectangular function Q(x) is unity if 0 < x ≤ 1 and vanishes otherwise 30 . 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 The bottom-left panel shows the semi-coherent detection statistic S as a function of f andḟ , maximized over the sky-grid. The value of S is represented by the colorbar. A further maximization over f is shown to the right; a further maximization overḟ is shown above. PSR J1620-4927, the darkest point near the bottom left of each figure, stands out clearly from the noise. In the same form, the right panel (b) presents the fully-coherent follow-up search results of the previous candidate (and every other candidate "dot") shown in the left panel. The quantity plotted is now the fully-coherent detection statistic P over the entire data set. As explained in the text, for each candidate this covers a region of parameter space which is four steps of the semi-coherent grid in each dimension.
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 an hierarchical scheme which is analogous to "zooming": successively swapping microscope objectives for ones of higher magnification, then re-centering 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) where for convenience we have normalized by the positive constant κ given by The computing cost to coherently follow up a single candidate is negligible in comparison to the cost of the previous semicoherent 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∆ f ) × (4∆ḟ ) × (4∆θ × 4∆θ) 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 , ∆ḟ , 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 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 which implies that the (complex) Fourier coefficients α are given by 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 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 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 spindown rateḟ selects the "narrowest and sharpest" overall pulse profile. In contrast, maximizing the statistic P = |α 1 | 2 of Equation (6) favors putting more power into the lowest harmonics.
Using H as the test statistic, a further stage of parameterspace 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 Hstatistic 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;Pletsch & Guillemot 2011). In the next section, the complete search scheme is illustrated with a detailed example.

EXAMPLE: RESULTS FOR PSR J1620-4927
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ḟ 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ḟ parameters.
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 parameterspace 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.
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 that 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 10 P FA ≈ −98. We have not tried to rigorously estimate the trials 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 21 , 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Ė = −4π 2 I fḟ , where I = 10 45 g cm 2 is the assumed neutron star moment of inertia, and the magnetic field strength at the neutron star surface B S ≈ 3.2 × 10 19 (−ḟ / f 3 s −1 ) 1/2 G and at the light cylinder B LC ≈ 2.94 × 10 8 (−ḟ f 3 s 5 ) 1/2 G, respectively. Figure 3 plots the newly discovered pulsars on a 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.

RESULTS FOR ALL NEW PULSARS
Two of the nine new systems, PSRs J1803-2149 and J2111+4606, are young energetic pulsars (Ė 6 × 10 35 erg s −1 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 11 G) of all LAT blind-search pulsars found to date. Also standing out, PSR J2139+4716 has the smallest spin-down power (Ė = 3 × 10 33 erg s −1 )  Table 1, all but one of the new pulsars have gamma-ray counterparts from the Fermi-LAT First Source Catalog (1FGL, Abdo et al. 2010a). The exception is PSR J2028+3332, which first appeared in the 2FGL catalog. In addition, these eight 1FGL unassociated sources have been classified as likely pulsar candidates based strictly on a comparison of their gamma-ray properties with previously detected LAT pulsars (Abdo et al. 2011a).
PSR J2111+4606 is associated with the source 0FGL J2110.8+4606 from the Fermi Bright Source List (Abdo et al. 2009b). Only ten objects from this catalog still remain to be identified.
The new system PSR J1620-4927 is associated with the unidentified AGILE source 1AGL J1624-4946 (Pittori et al. 2009). This AGILE source is however about ten times brighter, with a flux of (67 ± 13) × 10 −8 ph cm −2 s −1 . In addition, it has an error radius of 0. • 58 and encloses other 2FGL sources. The new gamma-ray pulsar could therefore only explain a fraction of the total flux detected by AGILE from 1AGL J1624-4946.
While not formally associated in the 2FGL catalog (Abdo et al. 2011b), PSR J2028+3332 lies within the 99% error contours for the EGRET source 3EG J2027+3429. However, the pulsar accounts for only about a fifth of the flux for the EGRET source. Another nearby LAT source, 2FGL J2025.1+3341, is associated with an AGN and has a peak LAT flux that can account for the remainder of the measured EGRET flux.
We have searched the TeVCat online catalog of TeV sources 31 for very high-energy counterparts to the newlydiscovered pulsars, but found no associations. This is not surprising, considering that most of the new gamma-ray pulsars have moderate spin-down power (∼ 10 34 erg s −1 ) and relatively large characteristic ages ( 100 kyr), while pulsars associated with TeV Pulsar Wind Nebulae (PWNe) tend to be young and energetic. PSRs J1803-2149 and J2111+4606 have properties that make associations with TeV PWNe plausible; however, no detections of TeV emission from these directions have been reported.
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 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 ( j w 2 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.) 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) singlepeaked 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 31 http://tevcat.uchicago.edu/ 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 singlepeaked 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 substructures 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 (cf. 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.

Spectral parameters
The spectral parameters for the new pulsars are obtained by fitting each phase-averaged spectrum with an exponentially cut-off 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 100 (in units of photons cm −2 s −1 ) and the energy flux G 100 (in units of erg cm −2 s −1 ) for events with energies between 100 MeV and 100 GeV.
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 is associated with a known supernova remnant, 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 spindown powerĖ for other gamma-ray pulsars with distance measures (cf. Abdo et al. 2010c). Based on this correlation "pseudo gamma-ray luminosities" L ps are derived as L ps ∼ 3.2 × 10 33 (Ė/10 34 erg s −1 ) 1/2 erg s −1 , where theĖ values are obtained from Table 2. Assuming a geometrical correction factor f Ω = 1 for the emission cone (Watters et al. 2009) for all gamma-ray pulsars, the relation L γ = 4π f Ω 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) 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? To examine whether the nine pulsars found with this new method could be detected with previous methods, we apply the same search method (Atwood et al. 2006;Ziegler et al. 2008) used to successfully discover the 26 previously-found blind-search LAT gamma-ray pulsars (Abdo et al. 2009a;Saz Parkinson et al. 2010. We select input data as done in previous searches; events are selected based on a fixed ROI and minimum-energy cut as described in Saz Parkinson et al. (2010). We use the same coherence window size (T = 2 19 s) as in the previous (and in the first stage of this paper's search ) search. No photon weights are computed or used. No skygridding is done in the first stage of the search: only the 2FGL-catalog sky position is used.
The previous blind search code recovers three of the nine pulsars (PSRs J1620-4927, J2028+3332 and J2111+4716) when the 2FGL sky locations are searched. Two of the remaining six pulsars are detected only if the correct known sky position is used (as opposed to the 2FGL-catalog position). If in addition the ROI and energy cuts are optimized (scanning different values) then three of the remaining four pulsars are detected. Finally, if the coherence window size is doubled (which dramatically raises the computational burden in a full blind search) then the last pulsar is found.
As compared to the previously-published blind-pulsar searches of LAT data, the methods used for this work incorporate several significant improvements in sensitivity and computational efficiency, as well as a longer data set, that explain why these new discoveries are made. First, the use of efficient parameter space gridding over both sky position and frequency derivative allows pulsars to be found that are much farther from the LAT catalog position than is possible with previous searches. In addition, using photon weights for both event selection and for the search computations ensures that the detection significance is near optimal with only a single trial over event selections. In contrast, methods that use a "cookie cutter" event selection must either search over two additional parameters (minimum energy and radius for the selection), or suffer a sensitivity penalty from potentially nonoptimal event selections. A key factor is that while weighting the photons provides only a modest sensitivity improvement over optimal cookie cutter selections, the improvement can be large when compared to non-optimal cuts. One might think that the additional systems found in this paper come about purely from the significantly-increased computer power that was available. This is not the case: the improved methods deserve almost all of the credit. The improved methods are so computationally efficient that had we searched only up to 64 Hz as was done in previous searches, we could have searched all 109 selected 2FGL sources over the initial spin-down range (|ḟ | ≤ 5×10 −13 Hz s −1 ) and would have found seven of the nine new pulsars using less than 5000 CPU hours. These are modest resources in comparison with those used in previous blind searches 32 . The remaining two pulsars would only have been found if the spin-down range were increased by a factor of twenty to |ḟ | ≤ 10 −11 Hz s −1 , increasing the required CPU time to about 100 000 hours. 6. RADIO AND X-RAY COUNTERPART 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 blindsearch 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 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 pulsarlike unassociated sources, all of these sources have 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 "trials 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): 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 3K 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 S ν ∝ ν α ). 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 −(θ/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 S 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 45minute GBT observations at 820 MHz, we detect the pulsar at a DM of 70.87 ± 0.2 pc cm −3 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 5 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 Green Bank Telescope 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 2 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 Pulsars and SNRs (Supernova Remnants) have the same origin, although the comparatively short lifetimes of SNRs means that the number of pulsar-SNR associations is quite small. However the rare associations are of high interest, as they constrain a number of pulsar and SNR parameters. A recent text-book example is the association of the gammaray pulsar J0007+7303 with the shell-type SNR CTA 1 (Abdo et al. 2008;Sun et al. 2011).
For the nine new gamma-ray pulsars no association with a known SNR listed in the most recent SNR-catalog (Green 2009) is found. However, SNRs are difficult to identify in case they are faint, confused or distorted by interaction with dense clouds. We have used large-scale and Galactic-plane radio continuum surveys to search for structures in the vicinity of the new pulsars, which may indicate an association. For the two pulsars PSRs J0106+4855 and J0622+3749 located above 10 • of Galactic latitude we used the MPIfRsurvey sampler 33 to extract maps from the 408 MHz and 1420 MHz surveys (Haslam et al. 1982;Reich 1982), which show faint sources and extended smooth diffuse emission, but no discrete object. The high-resolution interferometric NVSS maps (Condon et al. 1998) show numerous compact sources, but no extended features.
The remaining pulsars are located within 4 • of Galactic latitude covered by Galactic-plane surveys. The southern-sky pulsar PSR J1620-4927 is located towards the gradient of an extended emission complex as seen in southern-sky surveys (Reich et al. 2001;Jonas et al. 1998). The Southern Galactic Plane Survey (SGPS, Haverkorn et al. 2006) is insensitive to large scale emission, but shows no small scale structures within 0. • 5 of the pulsar. PSR J1746-3239 is located at b = −2. • 2 close to a emission ridge sticking out from the plane. The shell-type SNR G356.6-1.5 with a size of 20 ×15 (Gray 1994) is part of this ridge and about 0. • 6 away from the pulsar. PSR J1803-2149 is located just 2 away from the peak of a flat spectrum 6 Jy thermal source visible in all Galacticplane surveys. Quireza et al. (2006) list details for the HIIregion with a deconvolved Gaussian size of 2. 7 and a distance of 3.4 kpc. Also CO-emission was observed (Scoville et al. 1987). The pulsar is located at the periphery of this emission complex, but further studies are needed to investigate this region. PSR J2028+3332 is located at the southern boundary of the thermal Cygnus-X complex. The 11 cm survey (Reich et al. 1984) shows patchy structures in the field and an extragalactic 1.8 Jy source 40 distant from the pulsar. PSR J2030+4415 is seen towards the northern periphery of strong complex emission from the Cygnus-X region. Dedicated studies are needed to find any emission associated with the pulsar. PSRs J2111+4606 and J2139+4716 are both located in low emission areas, where the 21 cm and 11 cm Effelsberg Galactic-plane surveys show faint structures close to the noise level. 6.3. X-ray Subsequent to the pulsar discoveries, we searched for archival X-ray observations covering the new pulsars' sky locations. As listed in Table 7, we have found short Swift observations (3-10 ks exposure) for five of the pulsars. In addition there is a 6 ks-long XMM-Newton observation for PSR J1620-4927 and a 23 ks-long Suzaku observation for PSR J0106+4855. For the two remaining pulsars, PSRs J1746-3239 and J2028+3332, 10 ks-long Swift observations were carried out following the discoveries.
These X-ray data were analyzed, and no X-ray counterparts were found for any of the new gamma-ray pulsars.
In order to estimate upper limits on the X-ray flux for each new pulsar, a power-law spectrum with a photon index of 2 and a signal-to-noise of 3 is used. The values for the absorbing columns are estimated analogously to Marelli et al. (2011). In the 0.3-10 keV energy range, the derived upper limits on the X-ray flux of the pulsars are between 1 and 3×10 −13 erg cm −2 s −1 , except for PSR J1620-4927, where an upper-limit X-ray flux of 6×10 −14 erg cm −2 s −1 is obtained, and PSR J0106+4855, with an upper-limit X-ray flux of 2.3×10 −14 erg cm −2 s −1 . These results, listed in Table 7, are consistent with the gammato-X-ray flux-ratios for previously-found Fermi-LAT pulsars (Marelli et al. 2011). 7. CONCLUSION 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 pulsarlike 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 semicoherent, because coherent power computed using a window of 6 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, frequencyderivative 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 previouslyknown Fermi-LAT blind-search pulsar population (Abdo et al. 2009a;Saz Parkinson et al. 2010) by more than onethird, and brings 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 radioquiet 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 August 2011, we moved the computational burden of the search onto the volunteer distributed computing system Einstein@Home 34 (Abbott et al. 2009a,b;Knispel et al. 2010). This will provide significantly more computing power, and will allow a complete search of the parameterspace 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 leave 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 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 and 0970074.
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 operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.     NOTE. -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. Lps dps (GeV) (10 −8 photons cm −2 s −1 ) (10 −11 erg cm −2 s −1 ) (10 33 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.4 c 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 NOTE. -This table describes the spectral properties of each of the nine pulsars, modeling each spectrum as an exponentially cut-off power-law with photon indices Γ and cutoff energies Ec. 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 Lps and the pseudo distance dps are inferred from the apparent spin-down powerĖ and the energy flux G 100 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. a Photon flux measured above 100 MeV. b Energy flux measured above 100 MeV. c The 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. NOTE. -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. a Instrument-dependent sensitivity degradation factor. b The instrument records 800 MHz of bandwidth, but to account for a notch filter for RFI and the lower sensitivity near the band edges, we use an effective bandwidth of 700 MHz for the sensitivity calculations. NOTE. -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 Smin for each observation is computed at the observing frequency using Equation (14) and the parameters in Table  5, as described in the text. a Telescope pointing direction (not necessarily source position) b For these two observations, radio pulsations were detected at a flux density of 20 µJy (see text). c Observed 14 times for 1 hour each between this date and 2011-07-11. NOTE. -The sky locations of all nine pulsars have been searched for (non-pulsating) X-rays, using both archival and new data. No X-ray sources were found at the new pulsar locations, so flux upper limits and lower limits on the gamma-ray to X-ray luminosity ratio are reported. a Estimated analogously to Marelli et al. (2011). b Upper limit on the unabsorbed flux in the 0.3 − 10 keV energy range, using an absorbed powerlaw model with a photon index of 2 and a signal-to-noise of 3. c Lower limit on the Lγ/LX flux ratio.