EINSTEIN@HOME DISCOVERY OF 24 PULSARS IN THE PARKES MULTI-BEAM PULSAR SURVEY

The computing power for Einstein@Home is provided by members of the general public, donating idle compute cycles on their home and/or office PCs. These highly fluctuating computing resources are managed using BOINC software framework (Anderson et al. 2006) and a set of central servers administered by a handful of project developers and scientists.

BOINC is a framework to set up and manage a distributed volunteer computing project that requires minimal attention from the volunteers donating computing cycles. Members of the general public (volunteers) can register and attach their computers (hosts) to a BOINC project. For the volunteers, this only requires downloading¹³ and installing a small piece of software (BOINCManager) and selecting the distributed computing project of their choice. Currently, there are several dozens of BOINC-based distributed computing projects, covering a wide range of scientific disciplines.¹⁴

The BOINCManager collects information about the host and then coordinates the download of data files, analysis software, and processing instructions. Only computational tasks (work units) suitable for the given host are sent, and each work unit has a certain deadline (two weeks for Einstein@Home) after which a result has to be reported back to the project servers. Each work unit is sent to two independent hosts, the results are then compared based on predefined numerical metrics to ensure their scientific validity. If necessary, additional copies of the work unit are sent to other independent hosts until two agreeing results are found. A central database stores information relevant to the work units, volunteers, hosts, and internet discussion forums. Details about the internal BOINC processes controlling the work distribution as well as a more detailed discussion of the type of problems that can be solved by distributed volunteer computing projects are available in Allen et al. (2013).

Einstein@Home is one of the largest distributed computing projects. Since 2005, more than 330,000 volunteers have contributed to the project. On average, about 48,000 different volunteers donate computing time each week on roughly 150,000 different hosts. Currently, the sustained computing power is of the order of 1 PFlop s⁻¹. The radio pulsar search uses a varying fraction (between 20% and 50%) of the central processing unit (CPU) time available to the project. For the radio pulsar search, additional computing time is available on Nvidia and ATI/AMD graphics processing units (GPUs). Currently, approximately 10,000 hosts with Nvidia GPUs and 3600 hosts with ATI/AMD contribute to the project each week. For the PMPS analysis only executables for CPUs and Nvidia GPUs were available.

On average, 200 complete PMPS beams per day were analyzed by Einstein@Home between 2010 December and 2011 July, with three-day averages varying between 170 and 230 beams day⁻¹. The total computing time donated by the volunteers to analyze the data set described here is of the order of 17,000 CPU core years.

The PMPS data are publicly available and were copied from computer systems at the Jodrell Bank Observatory to the Albert Einstein Institute (AEI) in Hannover. For permanent storage at the AEI, data were copied to a Hierarchical Storage Management System, backed by a tape library. The filterbank data were also kept on spinning media to provide low-latency availability for pre- and post-processing purposes. In the following, we describe the preprocessing and dedispersion applied to each of the 41,364 observed beams in the Einstein@Home analysis.

The data were pre-processed on dedicated computers at the AEI Hannover. The survey data are provided in "filterbank format," i.e., radio frequency power spectra resolved in "filterbank channels," sampled at regular time intervals. In the first step, we convert the survey data in filterbank format from 1 bit dynamic range per sample into a file format with 8 bits per sample using tools from the SIGPROC software toolbox.¹⁵ This step is necessary for further processing by other software from the presto¹⁶ software toolbox (Ransom et al. 2002) and does not add dynamic range. The product of this first step is one filterbank file for each observed beam.

To mitigate the effect of RFI at later stages in the analysis, we used the presto software tool rfifind to obtain an individual RFI mask for each beam. We analyzed the 8 bit filterbank data in blocks of 2.048 s length in each filterbank channel and flagged persistent narrowband RFI and transient broadband RFI. In each beam, up to a few percent of the data were identified as RFI in this step. The next steps of the preprocessing replaced these blocks by constant values.

The data were down-sampled by a factor of two in the time domain, reducing the number of samples per time series to 2²² to save computational time. This was achieved by coadding neighboring time series bins, increasing the sampling time to 500 μs.

Free electrons in the interstellar medium delay the arrival time of the radio waves with a dependence on their frequency. If this effect is not corrected for, it smears radio pulses observed over a wide bandwidth, and severely reduces their detectability. The default method to mitigate this effect is to incoherently dedisperse the filterbank data, which is done by introducing frequency dependent delays to the different filterbank channels (Lorimer & Kramer 2005). The time delay relative to a reference frequency depends on the a priori unknown integrated electron column density along the line of sight, the dispersion measure (DM).

For the dedispersion with the presto tools, we used a set of 440 trial DMs up to 1876 pc cm⁻³, which exceeds the expected range of DM values in our Galaxy (Cordes & Lazio 2002). The set of trial values is chosen using the ddplan.py tool from presto and is shown in Table 1.

We used presto tools to dedisperse, barycenter, and down-sample the filterbank data. The result of the dedispersion step was a time series encoded with 32 bits per sample, and an additional header file for each DM trial value. A total of 440 dedispersed time series and corresponding header files was generated for each observed beam. Each of these down-sampled time series was 16.8 MB in size.

In the next step, we combined each dedispersed time series and its associated presto header file into a single file. The external header files are included as file headers containing relevant information about the time series data.

To save internet bandwidth for the data transfer to the Einstein@Home hosts, we encoded time series data with a dynamic range of 8 bits per sample. The 8 bit dynamic range is sufficient to encode the whole possible range of dedispersed time series samples. Given the original 1 bit sampling of the filterbank data and the number of filterbank channels (96), the maximum possible value of any de-dispersed time series sample is 96 ≈ 2^6.6 < 2⁸.

Each "compressed" time series had a total size of 4.2 MB. A bundle of four time series from a given beam formed the input data for a single Einstein@Home work unit. A single modern CPU core with an approximate computing power of ∼10 GFlops s⁻¹ can analyze each task containing 16.8 MB of data in ≈12 hr.

The ratio of data I/O to computing time is therefore of the order of 1 MB hr⁻¹, which means that a 30 MB s⁻¹ internet connection at a single download server can keep of the order of 10⁵ hosts continuously busy, assuming bandwidth is the limiting factor. Note that our search also runs on GPUs which can finish the same task about 10–20 times faster than a CPU; see Allen et al. (2013) for details. For GPUs the I/O to computing time ratio increases to up to 20 MB hr⁻¹.

A set of 400 preprocessed beams was kept ready for distribution at any time. This precautionary measure provided an ample time buffer for the case of technical difficulties (e.g., with the preprocessing machines) on the server side of the project.

In searching for possible pulsar signals hidden in instrumental noise, a signal model is required. Detection statistics are then designed to effectively identify that signal in detector noise. A full description of the signal model and detection statistics used in this paper are given in Allen et al. (2013). Here, we summarize the main points.

Our signal model describes the rotation phase Φ of the pulsar as seen in the radio telescope at time t, assuming that the pulsar is in a circular orbit:

$$\begin{equation} \Phi \left(t \right) = 2\pi f\left(t + \frac{a\sin \left(i\right)}{c} \sin \left(\Omega _{\rm orb} t + \psi \right)\right) + \Phi _0. \end{equation} \tag{ 1 }$$

Here, f is the intrinsic pulsar spin frequency and asin (i) is the projected orbital radius with inclination angle i. The orbital angular velocity Ω_orb is determined by the orbital period P_orb through Ω_orb = 2π/P_orb. The angle ψ denotes the initial orbital phase and Φ₀ is the initial signal phase.

To search for sinusoidal signals proportional to cos Φ(t), a large set of matched filters is applied to the instrumental output. Each matched filter is optimized for a particular waveform, and can be thought of as the "best possible search" for a signal at a particular point in the parameter space. The four parameters f, asin (i), Ω_orb,  and ψ are coordinates in this parameter space of possible signals. The detection statistics can be shown to be independent of the initial phase Φ₀; cf. Section 4.3 of Allen et al. (2013).

The matched filter at a particular point in parameter space will also respond to signals located "nearby," whose phase model is similar. Thus, in constructing a computationally efficient search, careful consideration must be given to the set of filters is chosen. The set of points in parameter space that is searched is called the template bank; Section 3.4 explains how it was selected.

At a particular point in signal parameter space, detection statistics are constructed from the average power $\mathcal {P}_n$ in the nth harmonic of the pulsar rotation phase Φ. In signal processing, the $\mathcal {P}_n$ are called "matched filter squared signal-to-noise ratios." This assumes Gaussian white noise with unit variance in the data stream. To search data for signals, the P_n are computed for first 16 harmonics of Φ (so n = 1, ..., 16) and then combined to form detection statistics.

The optimal way to combine the $\mathcal {P}_n$ into detection statistics depends upon the pulsar's pulse profile. For example, if the pulsar had a purely sinusoidal profile (at the rotation frequency) then the $\mathcal {P}_1$ would be the optimal detection statistic. However, since we are searching survey data for new pulsars, we do not know the pulse profile ab initio. Thus, we search five different detection statistics (often called harmonic sums) S₀, ..., S₄, constructed from summing different numbers of harmonics together:

$$\begin{equation} S_L \equiv \sum _{n=1}^{2^L} \mathcal {P}_n. \end{equation} \tag{ 2 }$$

The Lth harmonic sum is an optimal detection statistic (in the Neyman–Pearson sense, maximizing the detection rate for a given false-alarm rate) for a pulsar pulse profile that is a Dirac delta-function spike, truncated at the 2^Lth harmonic.

It is important to characterize the statistical properties of the S_L. In the absence of a pulsar, instrument noise can mimic a signal, resulting in large values of the detection statistics and a possible false alarm. If the instrument noise has Gaussian statistics, then in the absence of a pulsar it is easy to see that S_L behaves like a classical χ² random variable with 2^{L + 1} degrees of freedom. The false-alarm probability $p_{\rm FA}(S_L^*)$ is the probability that S_L exceeds a threshold value $S_L^*$ in the absence of a signal; for Gaussian noise this is given by an incomplete upper Γ-function.

We use the significance $\mathcal {S}$ to indicate the statistical significance of a signal candidate. Pulsar signals which are strong enough to be observable produce unusually large values of $\mathcal {P}_n$ and thus unusually large values of S_L. These in turn have low false-alarm probability. Hence, we define the significance of a candidate by

$$\begin{equation} \mathcal {S}\left(S_L\right) \equiv -\log _{10}\left(p_{\rm FA}\left(S_L\right)\right). \end{equation} \tag{ 3 }$$

For example, a candidate with significance $\mathcal {S} = 20$ has probability 10⁻²⁰ of occurring in random Gaussian noise.

The template grid or template bank is the set of points in the parameter space at which the detection statistic is evaluated. In an ideal world (unlimited computing power), the template bank would contain a very large number of signal templates: for any signal in the parameter space, there would be a template located nearby. No signal-to-noise ratio would be lost due to mismatch between the signal and template parameters.

Real template banks are designed to maximize detection probability at fixed computing cost. Substantial research work in the gravitational-wave detection community has shown how to construct (near-)optimal template banks (Owen 1996; Owen & Sathyaprakash 1999; Harry et al. 2009; Messenger et al. 2009; H. Fehrmann & H. J. Pletsch, in preparation). A real template bank is characterized by the its "worst-case" mismatch.

The mismatch m is defined (for a signal at one point in parameter space, and a filter template at a different point) as the fractional loss of detection statistic compared to a template placed at the signal location. In this paper, the same template bank is used for the detection statistics S₀, ..., S₄, however, the mismatch we discuss is only that of the "fundamental mode" $S_0 = \mathcal {P}_1$ detection statistic. The region of parameter space around a particular template, whose mismatch is less than some nominal value is called the template's coverage region. For small enough nominal values (m ≲ 0.01), the coverage region is an ellipsoid. However for the worst case mismatch used in this search (typically m = 0.2 or m = 0.3) the region is "banana shaped" as discussed in Section 3.5.

Our template bank and its construction method are similar to those described in Allen et al. (2013) and Knispel (2011). The four-dimensional template bank is obtained by a Cartesian product of a three-dimensional orbital template bank with a uniform frequency grid with spacing Δf = 1/3T. The extra factor of three is due to the mean padding used in our analysis; see Section 3.7 for details. The orbital template bank is not uniform. Since the region of constant mismatch varies as one moves over the parameter space, and is not ellipsoidal, the orbital template bank cannot be a regular lattice, constructed using standard methods such as (Owen & Sathyaprakash 1999). We use a combination of the random (Messenger et al. 2009) and stochastic (Harry et al. 2009) template-bank constructions. The orbital template bank was constructed using about 200k CPU hours on the Atlas cluster (Aulbert & Fehrmann 2009).

The resulting stochastic template bank has a nominal mismatch m₀ = 0.29 and coverage η = 90%. (The coverage is the fraction of points in parameter space whose mismatch m < m₀ from some template.) This efficient template bank contains 12,140 orbital templates; cf. about 60 acceleration trials in Eatough et al. (2013). The bank was tested using fake noise-free signals at grid points in orbital parameter space, with f = f_max, where f_max is the highest spin frequency for which the orbital template bank achieves the required nominal mismatch; see Section 3.5. As can be seen from Figure 1, the median mismatch m_med = 0.15 is much smaller than the nominal mismatch.

Zoom In Zoom Out Reset image size

To conduct a blind search for new pulsars, we must decide what region of parameter space to cover with a template bank. With unlimited computing power, we could search any parameter space, no matter how large. In practice, the finite computing power of Einstein@Home dictates that we only search some part of the parameter space.

Our choice is motivated first by astrophysical reasons: to target Einstein@Home to an interesting range of putative pulsar spin frequencies and orbital parameters. Second, for practical reasons, we decide to complete the Einstein@Home search on PMPS data within about a year, searching a previously unexplored part of the orbital parameter space. Thus, we constrain the search parameter space by setting a probabilistic limit on projected orbital radii, and by an upper limit on spin frequencies.

We estimated the available computing resources based on small-scale tests on a single computer and on the Atlas computing cluster at the AEI. The total Einstein@Home computing power was estimated from previous searches for radio pulsars.

The required number of orbital templates grows $\propto f_{\rm max}^3$. We constrained our orbital template bank to f_max = 130 Hz (P_min ≈ 7.7 ms).

Standard acceleration searches lose sensitivity where P_orb ≲ 10T. As discussed in Section 2, previous searches (Eatough et al. 2013) were sensitive for orbital periods P_orb ≲ 3 hr. Thus, we searched for orbital periods in the range 86 min ⩽ P_orb ⩽ 317 min. The upper boundary was chosen to provide seamless connection to the sensitivity ranges covered by previous searches. The lower boundary was dictated by the available computing power. Sensitivity to orbital periods as short as 86 min significantly increases the sensitivity to pulsars in compact binaries. We added a single template corresponding to an isolated pulsar to sustain sensitivity to isolated pulsars, or those in very wide orbits.

To provide full sensitivity along the complete orbit of any putative binary pulsar, the initial orbital phase ψ was not constrained, covering the range 0 ⩽ ψ < 2π.

We constrained the projected orbital radius by using a probabilistic bound on the orbital inclination angles based on the masses of putative pulsars and companions. From Kepler's third law, we define

$$\begin{equation} 0 \le a\sin (i) \le \alpha \frac{G^\frac{1}{3} \Omega _{\rm orb}^{-\frac{2}{3}}m_{\rm c,max}}{c(m_{\rm p,min} +m_{\rm c,max})^\frac{2}{3}}, \end{equation} \tag{ 4 }$$

where m_{c, max} is the maximum companion mass, m_{p, min} is the minimum pulsar mass, G is the gravitational constant, and c is the speed of light. The constant α measures the probabilistic bound on orbital inclination angles. For our search, we chose α = 0.5, m_{p, min} = 1.2 M_☉, and m_{c, max} = 1.6 M_☉. The scaling of the total number of templates with α, the pulsar and companion masses, and the orbital period range is non-trivial.

For different pulsar and companion masses, Equation (4) defines an upper limit on projected orbital radii as a function of the orbital angular velocity. For large companion masses, only a fraction of all possible orbital inclinations fulfill Equation (4). In other words, for smaller companion masses, the condition (4) is fulfilled for any i, and all of such binary pulsar systems are detectable by Einstein@Home.

Figure 2 shows the slice of orbital parameter space defined by these conditions. This choice guarantees that circular orbits at mass ratios for almost all pulsar-white dwarf binaries and many double neutron star systems lie inside our search space.

Zoom In Zoom Out Reset image size

Because we use the same template bank for all spin frequencies, it "overcovers" the parameter space for frequencies f < f_max. The required density of templates is highest at f_max, because the detection statistic depends upon the difference in phase, which varies most rapidly at the highest frequency. For frequencies higher than f_max, it "undercovers" the parameter space. In Figure 2, we show cuts through surfaces of constant mismatch around a single template with fixed orbital parameters and varying spin frequency. At lower frequencies, the region covered by the templates is "banana shaped" and extends well beyond the parameter space boundaries defined above. This is why our search detected the relativistic pulsar J1906+0746 (see Section 3.8) although its orbital parameters lie outside the parameter space covered by our template bank.

Our search loses sensitivity to the higher harmonics of fast-spinning MSPs, and searching for these higher harmonics is prohibited by the high computing costs. Increasing the value of f_max to (say) 1 kHz, is computationally unfeasible and would require ≈500 times more computing resources than our search!

As discussed in Section 3.5, standard acceleration searches (Ransom et al. 2002) lose sensitivity for binary pulsars with P_orb ≲ 10T. To quantify this effect for the PMPS data, we numerically computed the mismatch of an acceleration search and that of the Einstein@Home search at different orbital parameter space points and at the highest spin frequency f_max = 130 Hz. The orbital parameter space was covered by a cubic grid of 10,920 points in Ω_orb, asin (i), and ψ.

For each of the orbital parameter space points, we numerically simulated an acceleration search and the Einstein@Home search for a signal at this point. The simulated acceleration search followed the standard setup (Lorimer & Kramer 2005). We set up a grid of accelerations a₁ in a range of ±500 m s⁻² and in steps of δa₁ = 0.26 m s⁻². The acceleration range was chosen to agree with Eatough (2009), and the acceleration step size was chosen (very conservatively) requiring that the signal does not drift by more than half a Fourier bin as described in Section 6.2.1 of Lorimer & Kramer (2005). Further, we simulated the same frequency resolution Δf = 1/(3T) as for the Einstein@Home search. Then, we generated noise-free sine-wave signals given by Equation (1), and computed the detection statistic $\mathcal {P}_1 = S_0$, but using a template phase model Φ(t) which is a quadratic function of t, corresponding to the signal from a single-harmonic isolated pulsar moving with constant acceleration.

For the simulated Einstein@Home search, the detection statistic $\mathcal {P}_1=S_0$ was computed using a template phase model from Equation (1), then the fractional loss of detection statistic (mismatch m) was computed for both the acceleration and the Einstein@Home search.

Figure 3 shows the results of our comparison, where we have averaged the mismatch over the orbital phase ψ. For the simulated Einstein@Home search (left panel), the averaged mismatch $\overline{m}$ is close to the target value of nominal mismatch (m₀ = 0.3) in the entire parameter space. Because of the underlying random template bank, there are small parts with slightly higher and lower mismatches. Close to asin (i) = 0, the mismatch is considerably smaller, since this parameter space region is covered by the single template for isolated pulsars.

Zoom In Zoom Out Reset image size

For the simulated acceleration search (right panel of Figure 3), the mismatch reaches unacceptably high values $\overline{m} \gtrsim 0.5$ over a large fraction of the parameter space. Only in a small part does the acceleration search achieve mismatches comparable to the our method. Our results clearly show the improvement in the detection of radio pulsars in compact binary orbits. Regions of signal parameter space that were virtually inaccessible with acceleration searches are in reach of our method.

This section summarizes the "signal analysis" part of the Einstein@Home radio pulsar search pipeline, which runs on the volunteers' hosts and does the bulk of the computing work. The code is distributed under the GPL 2.0 license and is available for CPUs and GPUs under Linux, Windows, and Mac OS X; a complete description may be found in Allen et al. (2013).

For each host, the input data are typically four dedispersed time series which are analyzed sequentially as described here. The search code computes the detection statistics S₀, ..., S₄ at each template-bank point in parameter space, and then returns back to the Einstein@Home server a list of "top candidates": the points in parameter space where the detection statistic was largest. The analysis consists of a data preparation step, a loop over orbital templates, and an output/candidate reduction step.

In the data preparation step, the input time series is uncompressed and converted into single-precision floating-point format. It is whitened in the frequency domain using a running-median average spectrum. RFI is replaced by computer-generated Gaussian noise using a list of known contaminated frequency bands in the fluctuation power spectra of the zero-DM dedispersed PMPS data (Faulkner et al. 2004); this "zapping" step uses same list globally for the analysis of all PMPS beams. After this whitening and cleaning, the data are returned to the time domain.

Then the search code begins to iterate through the orbital templates. For each orbital template, the detection statistics S₀, ..., S₄. of Equation (2) are computed on the full frequency grid with spacing Δf = 1/3T. This is done by resampling the data in the time domain to remove the effects of the orbital modulation. The data are then mean-padded¹⁷ to length 3T and Fourier transformed using a fast Fourier transform (FFT); the detection statistic $\mathcal {P}_n$ is the squared modulus of the Fourier amplitude of the resampled time series in the nth frequency bin.

Inside the loop, the client search code maintains five different lists of candidates corresponding to the detection statistics S₀, ..., S₄. Each contains the 100 candidates with the largest values of S_i having distinct values of fundamental frequency f. The code checkpoints after each loop iteration; if needed, this permits it to be restarted with little loss of compute time.

When the loop over orbital templates is complete, the statistical significance $\mathcal {S}$ is computed for all 500 candidates, and the 100 with the highest $\mathcal {S}$ are returned to the Einstein@Home server in a single result file. Based on false-alarm statistics, one can show that selecting the 100 most significant candidates results in digging down well into the noise-dominated regime.

The result file is uploaded to the central Einstein@Home servers and validated in a two-step process, described in Allen et al. (2013). In order to be accepted, it has to agree with the result calculated by another volunteer's computer with a fractional accuracy of about one part in 10⁵.

After the upload of all 440 result files for each beam, one has to filter out the most promising candidates from all 44,000 candidates in each beam. The majority of candidates in any given beam are caused by random noise fluctuations and have low significance values $\mathcal {S}$. Thus, thresholding on $\mathcal {S}$ is a possible step in reducing the number of candidates. However, even in the presence of a highly significant signal, correlations can cause the signal to show up at multiple DMs, frequencies, and/or orbital parameters. Thus, more sophisticated methods are required to reduce the number of candidates to follow up.

We used two methods to identify pulsar-like signals, employing frequency coordinates tailored to the binary parameter space. One method is based on producing overview plots that allow a quick and easy identification of pulsars and first estimates of their orbital parameters. The second method employs automated filtering algorithms to identify promising candidates in the binary parameter space. Figure 4 shows a flow diagram comparing both methods, which we describe in the following.

Zoom In Zoom Out Reset image size

Our first method for the identification of promising pulsar candidates in the Einstein@Home result files uses custom-made overview plots. These visualize the complete set of candidates for any given PMPS observation. We identified pulsars by characteristic patterns in these plots. Promising candidates are followed up with tools from the presto software suite.

The set of overview plots is automatically produced for visual inspection when all valid result files for a given beam are available on the Einstein@Home servers. The plots show the significance for all 44,000 candidates as a function of a the spin frequency at the detector ν₁, the associated spin frequency derivative ν₂, and combinations of the orbital template parameters. The coordinates ν₁ and ν₂ resolve some of the correlations of the detection statistic $\mathcal {P}_n$ in the orbital parameters. The ν₁ and ν₂ are the linear and quadratic coefficients in a polynomial expansion of the phase model (1) in t. Their derivation may be found in the Appendix.

An example of the five different plots in our post-processing is given in Figure 5, which shows the highly significant detection of the binary pulsar J1906+0746 in the Einstein@Home results.

Zoom In Zoom Out Reset image size

The number of candidates identified from these overview plots is relatively small (a few per beam at most, none in the majority of beams, and on average 1 candidate in ∼50 beams). Tools from the presto software suite are used to fold the full-resolution filterbank data at the candidate spin period and DM identified from the Einstein@Home results and to optimize the parameters. In total, we inspected of the order of 1000 of these folded time series plots by eye and used them to judge the broadband nature and temporal continuity of the signal. The ATNF catalog¹⁸ (Manchester et al. 2005) and Web sites listing known, but as yet unpublished, pulsars¹⁹ are checked to ensure that the candidate is not a detection of a known pulsar. About half of our new discoveries were found by the visual inspection of overview plots.

The visual inspection of overview plots is augmented by an automated post-processing stage. This sifts through all candidates in a given beam and identifies the most promising ones for a follow-up via folding of the filterbank data. With infinite computing and man power, all candidates of a given beam would be used for folding the filterbank data. In practice, this is neither feasible nor necessary. Folding the 44,000 candidates for any given beam requires an unfeasibly large amount of computing time and human resources for further inspection. Many of the candidates will be caused either from random noise or by a single pulsar (or RFI) signal. For example, following up candidates at harmonically related frequencies and nearby DM trials does not yield additional information (cf. Figure 5).

The first goal of the automated post-processing is to reduce the number of candidates from the complete results files. The second goal is to find those pulsars whose identification with the overview plots is difficult. Weak signals at higher spin frequencies may only be registered at a few trial DMs above the noise level, making it difficult to detect them visually.

The automated post-processing consists of three steps. In the first step, we reduce the number of harmonically related candidates. We start at the candidate with the highest value of $\mathcal {S}$. Candidates associated with sub-harmonic and harmonic frequencies of this candidate are flagged as possibly related as follows. We conservatively assume a frequency band around the candidate's value of ν₁ with a width δν₁ given by the maximum Doppler modulation of all templates in the orbital template bank. The width of the frequency band is obtained from maximizing the frequency modulation amplitude Ω_orbasin (i)/c over the entire orbital template bank and computing the Doppler range at frequency ν₁ from

$$\begin{equation} \delta \nu _1 = \nu _1 \max (\Omega _{\rm orb}a\sin (i)/c). \end{equation} \tag{ 5 }$$

After this step, only candidates in the Doppler range around the fundamental frequency are retained for the next step.

These are further winnowed down in the second step. Real pulsar candidates should produce a peak in significance as a function of DM, separated from the noise-dominated background. To take this property into account, we approximate the significance as a function of DM (Cordes & McLaughlin 2003) by a parabola

$$\begin{equation} \mathcal {S}({\rm DM}) = \mathcal {S}_0 \left(1- \frac{{\rm DM} - {\rm DM}_0}{\Delta {\rm DM}}\right)^2. \end{equation} \tag{ 6 }$$

Here, $\mathcal {S}_0$ is the significance of the most significant candidate at central dispersion measure DM₀, and ΔDM is the "width" of the parabola, such that S(DM₀ ± ΔDM) = 0.

Now, we step over a fixed grid of 10 trial values of ΔDM between 0.6 pc cm⁻³ and 50 pc cm⁻³ and determine the best-fit parabola parameter ΔDM, i.e., the one with the smallest squared residuals. To make sure that the parabola approximation is accurate, we select candidates near the top of the peak with $\mathcal {S} \ge 0.6\mathcal {S}_0$. Out of all candidates within the Doppler frequency range and the best-fitting DM range we only keep the most significant candidate.

The procedure of removing harmonically related candidates and those at similar DMs described above is iterated with the next most significant candidate until only "independent" candidates are left.

The number of remaining candidates per beam was about 20 after this first reduction step. We then folded all reduced candidates using the presto software on the Atlas cluster at the AEI. Each candidate was folded at fixed DM with five different spin frequencies between ν₁ − 2δν and ν₁ + 2δν. Here, δν is the expected Doppler shift in spin frequency caused by orbital motion and is computed from the orbital parameters of the candidate. Folding at a conservatively wide parameter range accommodates possible offsets from the true orbital parameters. In this step, we folded of the order of 4,000,000 candidates for the whole PMPS search. Due to oversight, we folded the original data, not the ones resampled at the orbital parameters. We will rerun our post-processing and report results in a future paper.

The final step of our automated post-processing is done by an algorithm implemented in idl,²⁰ which selects pulsar-like candidates based on figures of merit derived from the prepfold plots and associated binary and ASCII files.

The algorithm first checks for disqualifying metrics in the ASCII file associated with the folded data output. Features that cause rejection include periods that fall within the range of known interference signals and a signal strength below a chosen threshold (5.0σ as computed by prepfold).

For output files that survive these tests, the software uses two different tests applied directly to the plots produced by prepfold. It applies them separately to two parts of each plot; see, e.g., Figure 7. The first test looks for signal strength and consistency in phase in the time versus phase plot (left-hand bottom subplot in Figure 7). A pulsar-like signal is a straight, vertical line when folded with the proper parameters. An average intensity is computed for each vertical line at every horizontal phase position by summing the pixel values in that column and dividing by the number of pixels summed. A pixel smoothing of 3, 9, and 21 phase bins is applied in the phase (horizontal) direction prior to taking each sum. This smoothing accounts for different possible pulse widths. The resulting average intensity values for the vertical columns are then compared to with an off-pulse average intensity value, and a different threshold test is applied for each smoothing case.

A second test looks at a contour plot showing signal strength as a function of trial period and period derivative (bottom right subplot in Figure 7). The expected signal is one that is localized with one significant peak. The algorithm counts and classifies the number of consolidated strong response regions (determined by the number of contiguous pixels that are saturated) into "large" and "small" areas, which have sizes of at least 8 and 250 pixels, respectively. These correspond to 0.03% and 0.86% of the total area of the bottom right subplot in Figure 7. The code retains the candidate as valid if 1 or 2 large areas and no more than 13 total large and small areas are counted, or if there are less than eight total large and small areas (regardless of the number of large areas). Violation of either of these conditions indicates that the plot is not likely to reflect a unique period and period derivative combination. Testing for both small and large areas was chosen to account for the fact that not all pulsar signals will be strong enough to yield just a single contour peak.

The thresholds described above were chosen by using a subset of the PALFA data collected with the WAPP back ends (Cordes et al. 2006) as a test case. A variety of area sizes and threshold numbers were tried for the tests to optimize the results (that is, until the number of candidates passing the test was reduced as much as possible without eliminating any known pulsars).

The idl code reduces the number of candidates by up to a factor of 100, leaving a number (∼100, 000) that we individually checked by eye.

Our two different pulsar identification methods have different selection effects, making them sensitive to different classes of pulsars. The automated processing has proven very useful in identifying signals which were missed in the overview plots, because they were weak and/or masked by RFI (cf. Table 2). But, since visual inspection of result plots and human judgment is the final step in both methods, human error is a possible source for the rejection of real pulsar signals. The simplified folding at five different spin frequencies as described above can lower the sensitivity to the most extreme compact binary pulsars in the automated candidate selection.

Table 2. Pulsars Discovered by the Einstein@Home Analysis of the Parkes Multi-beam Pulsar Survey Data

PSR	R.A.	Decl.	P	$\dot{P}$	P Epoch	DM	S1400	D	PMPS beam	$\mathcal {S}$	Discovery Method
	(J2000)	(J2000)	(s)	(10−15)	(MJD)	(pc cm−3)	(mJy)	(kpc)
J0811−38	08:11.7(5)	−38:57(7)	0.482594(2)	...	50824.5	336.2	0.3	6.2	0026_0051	15.6	V
J1227−6208♭	12:27.6(5)	−62:10(7)	0.034529685(8)	...	51034.1	363.2	0.8	8.4	0058_036D	17.9	A
J1305−66	13:05.6(5)	−66:39(7)	0.1972763(2)	...	51559.7	316.1	0.2	7.5	0109_0033	15.5	A
J1322−62	13:22.9(5)	−62:51(7)	1.044851(4)	...	50591.6	733.6	0.3	13.2	0001_0016	23.1	V
J1455−59	14:55.1(5)	−59:23(7)	0.1761912(2)	...	50841.7	498.0	1.6	7.0	0038_0182	14.0	V
J1601−50	16:01.4(5)	−50:23(7)	0.860777(4)	...	50993.6	59.0	0.4	3.6	0042_0039	29.1	A
J1619−42	16:19.1(5)	−42:02(7)	1.023152(4)	...	51975.6	172.0	0.6	3.7	0137_041B	35.4	V
J1626−44	16:27.0(5)	−44:22(7)	0.3083536(5)	...	51718.6	269.2	0.3	4.8	0125_077C	13.2	A
J1637−46	16:37.6(5)	−46:13(7)	0.493091(2)	...	50842.9	660.4	0.7	7.0	0039_0055	17.2	V
J1644−44	16:44.6(5)	−44:10(7)	0.1739106(2)	...	51030.2	535.1	0.4	6.2	0056_020B	14.1	V
J1644−46	16:44.1(5)	−46:26(7)	0.2509406(1)	...	50839.0	405.8	0.8	4.8	0035_0293	13.2	A
J1652−48♭	16:52.9(5)	−48:45(7)	0.0037851238(4)	...	51373.3	187.8	2.7	3.3	0085_0254	22.3	A
J1726−31♭	17:26.6(5)	−31:57(7)	0.12347018(9)	...	51026.4	264.4	0.4	4.1	0054_015A	15.9	A
J1748−3009♭	17:48:23.79(2)	−30:09:12.2(5)	0.009684273(2)	...	51495.1	420.2	1.4	5.0	0102_0059	18.0	A
J1750−2536♭	17:50:33.39(2)	−25:36:43(3)	0.034749053(8)	...	50593.8	178.4	0.4	3.2	0002_0089a	15.9	A
J1755−33	17:55.2(5)	−33:31(7)	0.959466(4)	...	52080.6	266.5	0.2	5.7	0141_0097b	21.2	V
J1804−28	18:04.8(5)	−28:07(7)	1.273011(9)	...	51973.7	203.5	0.4	4.2	0137_039B	13.2	A
J1811−1049†	18:11:17.07(8)	−10:49:03(4)	2.6238585620(3)	0.8(2)	55983.5	253.3	0.3	5.5	0149_0108	29.2	V
J1817−1938†	18:17:06.82(8)	−19:38.6(2)	2.0468376289(2)	0.36(9)	55991.8	519.6	0.1	8.6	0011_0323c	16.9	V
J1821−0331†	18:21:44.70(3)	−03:31:12.7(1)	0.90231562918(4)	2.53(2)	55980.9	171.5	0.2	4.3	0148_0197	28.3	V
J1838−01	18:38.5(5)	−01:01(7)	0.1832948(2)	...	51869.1	320.4	0.3	6.9	0132_0627	16.7	V
J1838−1849†	18:38:33.79(4)	−18:49:59(5)	0.48824200896(3)	0.04(1)	55991.9	169.9	0.4	4.5	0140_0064	31.7	V
J1840−0643♭†	18:40:09.44(5)	−06:43:47(1)	0.0355778755(3)	0.2202(7)	55930.0	500.0	1.2	6.8	0060_0206d	18.2	V
J1858−0736	18:58:44.3(7)	−07:37(7)	0.551058591(2)	5.06(7)	56108.5	194.0	0.3	5.0	0143_0051	16.7	A

Notes. Sources with a fully determined timing solution are marked with "†" and pulsars in binary systems by "^♭." For these sources, the values in parentheses are the 1σ errors as reported by tempo. For pulsars without a coherent timing solution, right ascension (R.A.) and declination (decl.) are the beam center coordinates or weighted averages in case of detections in multiple beams. Their position error is estimated from the PMPS beam size. Their spin periods P and dispersion measure DM are the values in the discovery observations. S₁₄₀₀ denotes flux densities and D the estimated distance. The discovery PMPS beam and the Einstein@Home significance $\mathcal {S}$ are given for reproducibility reasons. The last column shows the discovery method: "V" for the visual inspection, or "A" for the automated method. ^aJ1750−2536 was independently detected in the PMPS beam 0002_0096. ^bJ1755−33 was independently detected in the PMPS beams 0136_0268 and 0118_021A. ^cJ1817−1938 was independently detected in the PMPS beam 0043_0014. ^dJ1840−0643 was independently detected in the PMPS beam 0087_0026.

Download table as:  ASCIITypeset image

This pulsar in a binary system with a 6.7 day orbit was independently discovered by Mickaliger et al. (2012) and will be further characterized by D. Thornton et al. (in preparation). Although its significance $\mathcal {S} = 17.9$ in the Einstein@Home pipeline is relatively high, it is rather inconspicuous in the overview plots because of strong remnant RFI. The automated post-processing successfully identified this pulsar.

The discovery observation of this isolated 1045 ms pulsar at DM ≈734 pc cm⁻³ exhibits signs of intermittency. Follow-up observations confirmed this trend: we observed the pulsar three times at MJD 55592.7 (T ≈ 2100 s), MJD 55615.7 (T ≈ 590 s), and MJD 55648.6 (T ≈ 2100 s), respectively, with the Parkes telescope. We only detected the pulsar in the second (shortest) observation. The non-detections in the two other observations correspond to upper limits on the flux density of S ≲ 0.07 mJy, as computed from, e.g., Appendix 2.6 in Lorimer & Kramer (2005). Here and for all following upper limits we assumed a signal-to-noise threshold of 8σ, the observed pulsar duty cycle, sky temperature at the pulsar sky position, and system parameters as given in, e.g., Manchester et al. (2001).

The pulsar apparently shows long-term nulling or intermittent emission behavior (Kramer et al. 2006a; Lorimer et al. 2012). Like the other intermittent pulsars presented here, their large inferred distances mean that interstellar scintillation is unlikely to be the cause of the observed intensity variations; at this observing frequency scintles of width of only a few MHz are expected and therefore average out over the 288 MHz band. Further observations are required to quantify this effect for PSR J1322−63.

This isolated 176 ms pulsar showed a strongly pronounced scattering tail in its pulse profile in the discovery observation. Its DM ≈498 pc cm⁻³. Further measurements of the pulse shape at multiple radio frequencies could be used to study the interstellar medium, by determining exact values of the pulse-broadening timescales τ_d. Assuming a delta-function pulse shape, we determine τ_d = 0.08(4) s from the discovery observation. This is in very good agreement with the measurements in Bhat et al. (2004). More exact measurements of τ_d in turn can be used to improve Galactic electron density models (Cordes & Lazio 2002), and to update existing empirical relations between τ_d and the DM (Bhat et al. 2004).

This binary MSP has a spin period of 3.78 ms and a DM of ≈188 pc cm⁻³. This pulsar has the fifth highest value of DM/P = 49.7 pc cm⁻³ ms⁻¹ of all Galactic field pulsars. As discussed above, the ratio is a measure of the depth to which pulsar surveys probe our Galaxy for MSPs. The value for PSR J1652−48 is only surpassed by that for PSR J1903+0327, PSR J1900+0308, and PSR J1944+2236 found in the PALFA survey (Champion et al. 2008; Crawford et al. 2012), and by that for PSR J1747−4036, found through a radio search of an unidentified Fermi source (Kerr et al. 2012). The pulsar is being timed at the Parkes Observatory.

We confirmed PSR J1726−31 upon the first redetection attempt at Parkes telescope in 2012 February. In nine subsequent observations at Parkes, only one further detection of this pulsar at MJD 56111.5 has been made. The other non-detections were done with the following integrations times (upper limit on the flux density in parentheses): once with T ≈ 3600 s (S ≲ 0.08 mJy), four times with T ≈ 2100 s, once with each T ≈ 1800 s, T ≈ 1500 s, and T ≈ 1300 s (S ≲ 0.1 mJy in each case). The large distance of 4.1 kpc, inferred from DM ≈ 264 pc cm⁻³ means interstellar scintillation is unlikely to be the cause of the observed intensity variations.

A search of the Parkes archival observation logs has revealed that this pulsar was independently identified as a promising candidate in 1999, at the time of the original PMPS observations. Four confirmation attempts were made in the same year, but without success. Three of them have integration times T ≈ 2100 s and each limits the flux density to S ≲ 0.1 mJy, one has T ≈ 360 s and therefore S ≲ 0.25 mJy.

Including our observations with these earlier confirmation attempts, and the original survey observation, the pulsar has been observed for a total of ≈8.6 hr. Of this time the pulsar is visible for ≈1.9 hr, suggesting the pulsar is detectable about 20% of the time.

In addition to the intermittency, there is evidence for large period changes between observations, and one instance of significant line-of-sight acceleration (a ≈ 11 m s⁻²). This suggests that this pulsar is a member of a compact binary system. Additional timing observations to characterize this pulsar, its possible companion, and the orbital parameters are ongoing.

If this pulsar is in a compact binary system, it is strikingly similar to PSR J1744−3922 (Faulkner et al. 2004). PSR J1744−3922 is mildly recycled with a spin period of 172 ms (cf. 123 ms for PSR J1726−31) and resides in a compact binary system (P_orb ≈ 4.6 hr). Like our discovery, PSR J1744−3922 also exhibits nulling and is visible for approximately 25% of the time at 1.4 GHz.

PSR J1744−3922, and thus PSR J1726−31, might be the members of a new class of binary pulsars as proposed by Breton et al. (2007): these binary pulsars have long spin periods, large magnetic fields (∼10^{10 − 11} G), low-mass companions, and low orbital eccentricities. Their evolutionary history is not well understood and known formation channels fail to explain all of their properties.

A different possible explanation for the observed intensity variations are eclipses in a compact binary systems: "black widow" systems like PSR B1957+20 (Fruchter et al. 1988). This seems unlikely, though, since the spin period of PSR J1726−31 is significantly larger than those of other known eclipsing binary pulsars (Freire 2005).

This pulsar is a 9.7 ms pulsar in a 2.93 day binary orbit. With DM ≈ 420 pc cm⁻³ it has by far the highest DM of all MSPs known to date. Its value of DM/P = 43.3 pc cm⁻³ ms⁻¹ is the seventh highest value published (Manchester et al. 2005; Crawford et al. 2012). The pulsar is currently being observed at Jodrell Bank on a regular basis to determine a full timing solution, which we will publish in a second paper on our search. Table 3 shows first measurements of its orbital parameters. The mass function f ≈ 2.87 × 10⁻⁴ M_☉ of this system indicates a minimum (median) companion mass of 0.09 M_☉ (0.10 M_☉) for a pulsar mass of 1.4 M_☉.

This pulsar has a spin period of 34.7 ms and was discovered independently in two adjacent PMPS beams. It has a DM of ≈178 pc cm⁻³ and is a member of a binary system with an orbital period of 17.1 days. It is currently being observed regularly at Jodrell Bank to improve the timing solution, which will be published in a follow-up paper. Table 3 shows first measurements of its orbital parameters. The mass function computed from the orbital parameters is f ≈ 0.029 M_☉, yielding a minimum (median) companion mass of 0.47 M_☉ (0.56 M_☉) for a pulsar of 1.4 M_☉.

Besides the "common" low-mass binary pulsars (LMBPs) exists the rather rare class of intermediate-mass binary pulsars (IMBPs). IMBPs differ from the LMBPs by longer spin periods, more massive companions, and larger orbital eccentricities. Further, their binary evolution channels seem to be significantly different (Camilo et al. 2001).

The combination of a long orbital period (17.1 days), relatively long spin period (34.7 ms), and an eccentricity of e = 3.92(4) × 10⁻⁴ makes PSR J1750−2536 most likely an IMBP. Plotting the orbital parameters of PSR J1750−2536 in the P_orb–e plane clearly places this pulsar outside the parameter space expected for LMBPs as predicted by a relation found by Phinney (1992; cf. Figure 4 of Camilo et al. 2001). Its low inferred distance from the Galactic plane |z| = 0.04 kpc is comparable to that of the other known IMBPs.

We will determine the spin-down of the pulsar with an ongoing timing campaign, from which we will in turn obtain the surface magnetic field of the pulsar. If PSR J1725−2536 is indeed an IMBP, its magnetic field B is expected to be relatively low at B = (5–90) × 10⁸ G.

This source was discovered by the Einstein@Home project in 2011 February in two independent PMPS beams and confirmed by a second observation in 2011 July. It has been independently discovered and published without a timing solution by Bates et al. (2012). We observed the pulsar regularly at Jodrell Bank and obtained a timing solution, which we report in Table 2. The discovery observation of this 2047 ms pulsar at DM ≈520 pc cm⁻³ exhibited signs of intermittency: it showed fading and re-appearing radio emission over the course of both 35 min discovery observations. We re-observed the pulsar with the Parkes telescope at MJD 55649.7 for T ≈ 1800 s and performed a gridding observation the Effelsberg 100 m telescope at MJD 55725.0. The pulsar was not detectable in the Parkes observation (thus S ≲ 0.09 mJy), but we found it permanently emitting in an 11.5 min observation with the Effelsberg telescope.

This 902 ms pulsar with DM ≈172 pc cm⁻³ exhibited signs of intermittency in its discovery observations. A follow-up observation on MJD 55790.6 using the Parkes telescope did not convincingly confirm the pulsar (T ≈ 1800 s, S ≲ 0.06 mJy), but a second observation on MJD 55808.2 at Jodrell Bank did. The first confirmation attempt and varying flux density in the Jodrell Bank observation confirm that this pulsar shows intermittent emission on long timescales. Table 2 shows its properties from a fully determined timing solution, obtained using observations with the Lovell telescope at Jodrell Bank.

This 35.6 ms pulsar was identified independently in two PMPS survey beams at DM ≈500 pc cm⁻³. The barycentric spin period observed in each one of these observations shows an increase that is inconsistent with pulsar spin-down. Ongoing timing observations at Jodrell Bank have revealed that this pulsar is in a binary system with an unusually large orbital period of 937 days (the fourth largest known) and with a small eccentricity, e ≈ 0. Table 3 shows measurements of its orbital parameters from our coherent timing solution. The mass function f ≈ 1.77 × 10⁻³ M_☉ for this system yields a minimum (median) companion mass of 0.16 M_☉ (0.19 M_☉) for a pulsar mass of 1.4 M_☉.

If the companion is a typical low-mass He white dwarf, as the mass function suggests, plotting this pulsar on the Corbet diagram of spin period versus orbital period (Corbet 1984) reveals that this pulsar is located in the same region as other long orbital period systems (P_orb > 200 days) which contain recycled pulsars with longer spin periods (Tauris 2011). Tauris & Savonije (1999) provide an explanation for the origin of such systems: the progenitor was likely a wide low-mass X-ray binary system. In such systems, there is only a short period of mass transfer because the donor star is highly evolved by the time it fills its Roche lobe.

If the companion is indeed a white dwarf this system is a potential target for tests of the strong equivalence principle via the Damour–Schaefer test (Damour & Schaefer 1991). Unfortunately, it can be shown that at least one of the fundamental criteria for this test is not fulfilled: the angular velocity of the relativistic advance of periastron, $\dot{\omega }$, must be appreciably larger than the angular velocity of the systems rotation around the Galactic center, Ω_Gal. For this system, and based on the DM-inferred distance, Ω_Gal ≈ 3.6 deg Myr⁻¹, and $\dot{\omega } \approx 3.0$ deg Myr⁻¹ (N. Wex 2012, private communication). If the system orientation is well understood, tests of SEP can be made via other methods, however, this also requires high timing precision (Freire et al. 2012a).