A Millisecond Pulsar Discovery in a Survey of Unidentified Fermi γ-Ray Sources with LOFAR

We have used the LOFAR High Band Antennas (HBAs; van Haarlem et al. 2013) of 21 of the 24 LOFAR Core stations¹⁰ to form 7 tied-array beams (Stappers et al. 2011). This observational setup has baselines up to 2.3 km and provides tied-array beams of  ∼35 in diameter (FWHM) at the central frequency, with 7 beams covering a total circular field-of-view (FoV) of about 10' in diameter (see Figure 1). With this setup we have observed 52 out of 1010 unidentified γ-ray sources from the 3FGL Fermi-LAT point-source catalog (Acero et al. 2015). These 52 sources were selected as they are visible to LOFAR (source elevation >30° during transit), located away from the Galactic plane ($| b| \gt 10^\circ ;$ where the sky temperature and scattering at 135 MHz are significantly lower), and because they have positional uncertainty regions less than 10' in diameter (i.e., fit the FoV of 7 tied-array beams). No cuts on the spectral parameters of the sources were performed. The observed sources and some of their parameters are listed in Table 1. The sample of Fermi sources searched here does not overlap with that of Bassa et al. (2017b).

Zoom In Zoom Out Reset image size

We employed a semi-coherent dedispersion scheme, aimed at mitigating the effects of dispersive smearing and implemented in cdmt (Bassa et al. 2017a). To allow coherent dedispersion, we have recorded complex voltage data for dual-polarization, Nyquist sampled subbands of 195.3125 kHz bandwidth (5.12 μs sampling). To maximize sensitivity and FoV we have used signals from 200 subbands in the 115−155 MHz frequency range (39.06 MHz bandwidth). Modest integration times of ${T}_{\mathrm{obs}}=20$ minutes were chosen to maintain sensitivity to accelerated signals from binary systems.

For each observation, the 200 frequency subbands were coherently dedispersed to 80 evenly spaced trial dispersion measures (DMs), ranging from 0.5 to 79.5 pc cm⁻³ (about twice the expected maximum Galactic DM for most of the surveyed sources), and channelized into a total of 1600 channels, using cdmt. The time and spectral resolution after channelization were 40.96 μs and 24.41 kHz, respectively. Around each coherent DM trial we made incoherent DM trials in steps of 0.002 pc cm⁻³. The two DM step sizes are chosen to limit the total (intra-channel and ΔDM) dispersive smearing compared to the true DM of the source to a maximum of 0.15 ms (see the top panel in Figure 2). Each dedispersed time series was searched for accelerated periodic signals in the frequency domain, and the 200 best pulsar candidates from each beam, according to a modified version of PRESTO's accel_sift.py sifting script (Ransom 2001), were folded and inspected by eye.

Zoom In Zoom Out Reset image size

Confirmation observations used all LOFAR Core stations (baselines up to 3.5 km), and thus have tied-array beams with a ∼3 times smaller area of ∼2' in diameter (FWHM) at the central frequency. Furthermore, the tied-array ring size (the offset of the center of the outer beams from the center of the pointing) is reduced to 175, such that the beams overlap slightly, and the position of a newly discovered pulsar can be refined by weighting the signal-to-noise ratios of detections in the different beams (see the colored filled circles in Figure 1 for an illustration). However, note that the ionosphere can shift beams by approximately an arcminute during periods of strong ionospheric turbulence. This can somewhat reduce the accuracy of this positional determination method.

Although the effects of dispersive smearing within a channel can be mitigated by the use of coherent dedispersion, the sensitivity of any pulsar survey at low radio frequencies is ultimately limited by scattering (approximately $\propto \,{\nu }_{\mathrm{obs}}^{-4}$), which results in an exponential broadening of the observed pulse shapes. We calculated the expected scattering times using the empirical fit for scattering as a function of DM made by Bhat et al. (2004), and compared this to dispersive smearing within channels (see the top panel in Figure 2). We have calculated the minimum detectable flux density our survey was sensitive to, using the modified radiometer equation for pulsars (Lorimer & Kramer 2012, Appendix 1.4), where we have used σ = 10 as the minimum signal-to-noise ratio for a convincing pulsar candidate (although candidates with a somewhat lower signal-to-noise ratio were also investigated), $\beta \,\approx$ 1.0 as the digitization correction factor (survey observations were processed with 8-bit integer bit depth), ${T}_{\mathrm{sys}}\,\approx$ 400 K as the temperature of the telescope and the sky at the observing frequency, and $G\,\approx$ 5.6 K Jy⁻¹ as the telescope's gain.¹¹ Sensitivity curves for a pulsar with a 1 and a 10 ms spin period and an intrinsic 10% duty cycle are shown in the bottom panel of Figure 2. Out to DMs of about 40 pc cm⁻³ we were sensitive to 2 ms pulsars brighter than ∼2 mJy, if the source was not eclipsed at the time of observation (many binary γ-ray MSPs are eclipsed for up to ∼50% of their orbit). This flux limit applies to observations at zenith; the sensitivity falls off approximately as ${\sin }^{-1.4}({\theta }_{{\rm{z}}})$, where ${\theta }_{{\rm{z}}}$ is the zenith angle (Noutsos et al. 2015).

Following the discovery and confirmation of the pulsar we started a timing campaign with LOFAR. Timing observations use all Core stations and the HBA bandwidth from 110 to 188 MHz, two times wider than possible in the survey observations. After initial dense and logarithmically spaced 10 minute observations spanning two weeks, the pulsar was observed once per month for 20 minutes. All observations are dedispersed and folded using the LOFAR Pulsar Pipeline (e.g., Kondratiev et al. 2016). Pulse times-of-arrival (TOAs) are extracted from 5 minute sub-integrations using tools from the PSRCHIVE¹² (Hotan et al. 2004) pulsar software package. We have used TEMPO2¹³ (Hobbs et al. 2006) to obtain an initial phase-connected timing solution spanning 0.74 years and fitted for position, spin frequency, and DM (see Figure 4). The efac/equad plug-in (Wang et al. 2015) was used to rescale the LOFAR TOA uncertainties, suggesting that a multiplication factor of 1.3 and an additional uncertainty of 0.8 μs (multiplied with and added to the initial uncertainty in quadrature) better reflect the expected Gaussian scatter of the residuals. Note that scattering can influence the measured DM and that there are thus likely systematic uncertainties on the DM that are larger than the nominal TEMPO2 error listed in Table 2. Also, the frequency dependence of the pulse profile might bias the measured DM value.

Zoom In Zoom Out Reset image size

Table 2.  Parameters for PSR J1552+5437

Parameter	Value
Timing Parameters (Radio and γ-Ray)
R.A. (J2000)	15h52m$53\buildrel{\rm{s}}\over{.}$33117(17)
Decl. (J2000)	+54°37'057866(14)
Spin frequency (Hz)	411.88053142429(10)
Frequency derivative (Hz s−1)	−4.746(17) × 10−16
Dispersion measure (pc cm−3)	22.9000(5)
Span of timing data (MJD)	54871.7–57698.5
Epoch of timing solution (MJD)	56285
Number of TOAs	88
rms timing residual (μs)	10.1
Reduced ${\chi }^{2}$ value	1.1
Clock correction procedure	TT(BIPM2011)
Solar system ephemeris model	DE421
Radio Flux Densities
Flux density at 150 MHz (mJy)	3.8 ± 1.9
Flux density at 820 MHz (μJy)	<17
Flux density at 1.4 GHz (μJy)	<20
Derived Parameters
Spin period (ms)	2.4279
Spectral index	$\lt -2.8\pm 0.4$
Galactic longitude (°)	85.6
Galactic latitude (°)	47.2
DM-derived distancea (kpc)	1.2, 2.6
Spin-down luminositya,b (erg s−1)	(8.9, 9.4) × 1033
Surface magnetic fielda,b(G)	(9.0, 9.2) × 107
Characteristic agea,b(years)	(1.2, 1.1) × 1010
γ-Ray Parameters
γ-ray-radio profile lag (ϕ)	0.042 ±  0.004 ± 0.1
γ-ray peak separation (ϕ)	0.53 ± 0.01
γ-ray photon index	1.4 ± 0.3
γ-ray cutoff energy (GeV)	3.7 ± 1.6
Photon flux (cm−2 s−1)	(2.5 ± 0.9) × 10−9
Energy flux (erg cm−2 s−1)	(2.7 ± 0.4) × 10−12
Luminositya(1032 erg s−1)	(4.7 ± 0.7), (22 ± 3.2)
Efficiencya(%)	6.1, 28.6

Notes.

^aBased on the NE2001 (Cordes & Lazio 2002) and the YMW16 (Yao et al. 2017) models, respectively. ^bUpper limit: corrected for the acceleration due to the kinematics of the Galaxy, but not for the Shklovskii effect.

Download table as:  ASCIITypeset image

The pulsar's flux density was measured in all timing observations by calibrating the observations using an improved Hamaker beam model (Hamaker 2006) and comparing the on-pulse with the off-pulse window (full details of LOFAR MSP flux calibration are described by Kondratiev et al. 2016). These measurements lead to a mean flux density for 19 observations at 150 MHz of 3.8 ± 1.9 mJy (50% uncertainty), but the observed flux density can vary by a factor ∼2 from observation to observation—possibly because of refractive scintillation, though RFI and ionospheric beam jitter can also influence this. A search for the Faraday rotation measure toward the pulsar using PSRCHIVE's rmfit routine did not converge for any of the LOFAR observations, likely because the pulsar shows little or no polarization beyond the detection limit.

PSR J1552+5437 was observed at L-band for $9\times 30\,\mathrm{minutes}$ and 4 × 1 hr with the Lovell (400 MHz bandwidth at 1532 MHz center frequency) and Nançay (NRT; 512 MHz at 1486 MHz) radio telescopes, but not detected, limiting the flux density to less than 13 μJy, when these data sets are co-added (under the assumption that diffractive scintillation is averaged out; here, the expected scintillation bandwidth is only ∼10 MHz at 1400 MHz). A non-detection of the pulsar at 820 MHz (200 MHz bandwidth) in a 1.5 hr observation using the Robert C. Byrd Green Bank Telescope sets an upper limit to the pulsar's flux density at 820 MHz of 17 μJy. These upper limits are calculated using the modified radiometer equation for pulsars (Lorimer & Kramer 2012, Appendix 1.4), assuming that PSR J1552+5437 would have been detected in those bands if it had a signal-to-noise ratio of at least 5 and a similar pulse width to our LOFAR detections. Based on the same assumptions, the upper limit for NRT is confirmed to be 20 μJy after flux calibration of the 4 × 1 hr of observation with a pulsed noise diode and a calibration source. The detections at 150 MHz and the upper limits at 820 and 1400 MHz constrain the radio power-law spectral index of PSR J1552+5437 to be $\alpha \lt -$ 3.2 ± 0.4, where we assumed that the radiometer equation has a 50% uncertainty. A more conservative upper limit, however, takes into account the potential overestimate of LOFAR pulsar fluxes by a factor ∼2 that was noted by Frail et al. (2016). In that case the LOFAR 150 MHz flux density would be 1.9 ± 0.9 mJy, and the upper limit on the spectral index $\alpha \lt -$ 2.8 ± 0.4. Future observations at 350 MHz will also be useful for mapping the spectrum.

We have also observed the pulsar with LOFAR's Low Band Antennas (LBAs) for 1 hr at $30\mbox{--}90\,\mathrm{MHz}$ on MJD 57496, but were unable to detect the pulsar by folding the data with the best-fit parameters derived from the timing analysis. This is unsurprising given the faintness of the source in the LOFAR HBA and the increased system temperature T_sys in the LBA. Also, the scattering tail that is already visible in the radio profile at 150 MHz (Figure 3) will be $\sim 10\times$ larger in the LBA range and would smear out the pulsations. Only three (very bright, and unscattered) MSPs have so far been detected using the LOFAR LBAs (Kondratiev et al. 2016).

We downloaded the Fermi-LAT Pass 8 photons of the SOURCE class from 2008 August 4 (the start of the mission) to 2016 October 14, within $20^\circ$ of the best position derived from radio timing. Using the Fermi Science Tools, we selected the photons in the energy range 0.1–100 GeV using the recommended cuts. We performed a binned maximum likelihood gtlike analysis on the photons in the 20° × 20° square centered on the timing position, leaving only the spectral parameters of the sources within the inner 5°  radius free. Our source model was based on the 3FGL catalog and as models for the Galactic diffuse emission and isotropic diffuse background we used the gll_iem_v06.fits (Acero et al. 2016) and iso_P8R2_SOURCE_V6_v06.txt templates¹⁴ , respectively. 3FGL J1553.1+5437 moved to the pulsar's timing position is detected with a test statistic TS value of 205 (about 14σ, while the source had a ∼ 8.5σ significance previously) using an exponentially cutoff power-law model to describe its spectrum. The exponentially cutoff model is preferred over a simpler power law as TS_cut ≡ 2Δ log(likelihood) = 15 > 9 (likelihood ratio test, following Abdo et al. 2013), and the best-fit parameters are listed in Table 2.

Based on the spectral analysis, all the events in the region around the source were assigned a probability of originating from 3FGL J1553.1+5437 using gtsrcprob (Kerr 2011). Selecting only those events with a probability $\gt 20$% resulted in 350 photons. Pulsar rotational phases ${\phi }_{i}(t)$ were computed based on the radio timing solution using TEMPO2¹⁵ (Hobbs et al. 2006) with the fermi plug-in (Ray et al. 2011). Folding the γ-ray photons over the range where the radio timing solution was valid did not result in a significant pulse profile, and we thus performed a brute-force search over the pulsar's spin frequency f and spin-frequency derivative $\dot{f}$ to find a coherent solution over the 7.5 years of Fermi data (neglecting higher order effects in this search is feasible because the MSP is likely isolated).

In the brute-force search, the barycentered phases were updated using the Taylor series

$$\begin{eqnarray}&&{\phi }_{i}(t)={\phi }_{i,0}+f({t}_{i}-{t}_{0})+\displaystyle \frac{1}{2}\dot{f}{({t}_{i}-{t}_{0})}^{2}\end{eqnarray} \tag{ 1 }$$

for 100 × 100 values of f and $\dot{f}$ within two times the error range of the radio timing solution. The H-test (de Jager et al. 1989) of the folded pulse profile was calculated for each trial. With the f and $\dot{f}$ that maximized H to 70, it was possible to significantly fold all 350 Fermi photons, which confirms the link between PSR J1552+5437 and 3FGL J1553.1+5437.

To lift the degeneracy between astrometric and rotational parameters in the timing solution we included the γ-ray data in our timing analysis. We used an unbinned maximum likelihood method to extract 8 topocentered TOAs with at least a 3σ detection from the 350 Fermi photons (Ray et al. 2011). More sophisticated and sensitive unbinned methods for extracting γ-ray TOAs have been developed in recent years (e.g., Kerr et al. 2015; Pletsch & Clark 2015), but using the method described above suffices for our present purposes.

The results of joint radio and γ-ray timing are listed in Table 2, the γ-ray profile folded with the final timing solution is depicted in Figure 3, and the timing residuals as a function of time are shown in Figure 4. The timing position is not at the center of the three beams with the best detections in the confirmation observations (see Figure 1), while the timing position based on the radio data alone agrees with the full timing solution to within a few hundredths of an arcsecond. Possibly, the ionosphere has caused the beams to shift by ∼1' in the confirmation observation. We also note that the schematic shown in Figure 1 is only a rough approximation of the true beam shapes.

The observed spin period derivative of 2.80 × 10⁻²¹ s s⁻¹ is not the intrinsic value, as it has to be corrected for the non-zero proper motion of the pulsar, the Shklovskii effect (Shklovskii 1970), and for movement due to the kinematics of the Galaxy (e.g., Nice & Taylor 1995). The Galactic contribution is (−4.36, −6.20) $\times \,{10}^{-22}$ s s⁻¹ for a (1.2, 2.6) kpc distance in the line-of-sight. This is the sum of the differential Galactic rotation and the k_z term. Adding this correction leads to a spin-frequency derivative of (3.24, 3.42) × 10⁻²¹ s s⁻¹. With the current data, it was not possible to significantly fit for the proper motion of the pulsar. However, the uncertainty on the fit values limits the proper motion to <36.8 mas yr⁻¹ (3σ), corresponding to a Shklovskii correction to the spin period derivative of <(9.6, 20.8) $\times \,{10}^{-21}$ s s⁻¹.The inferred surface magnetic field strength based on the observed spin period and spin period derivative, corrected for Galactic acceleration, is with (9.0, 9.2) × 10⁷ G already one of the lowest pulsar magnetic fields measured to date, and will become slightly lower after correcting $\dot{P}$ with an extended timing baseline.

Finally, we consider the offsets between the radio and γ-ray pulse peaks. As can be seen in Figure 3, both the radio and the γ-ray profile show a main pulse and a subpulse offset by about half a rotational phase. The γ-ray profile does not show any additional features when the number of phase bins is increased. We have set the rotational phase $\phi =0$ at the onset of the main pulse of the LOFAR radio profile. To quantify the peak separations, we fitted a Gaussian profile to both radio pulses. The peak around phase 0 in the probability-weighted γ-ray profile was fitted using two Lorentzian profiles, and the other pulse with one Lorentzian profile, on top of the background. The maximum of the radio profile is at ${\phi }_{{\rm{r}}}=0.063\pm 0.002$ (where the rotational phase is defined between 0.0 and 1.0, and errors are statistical), the radio subpulse peaks at $\phi =0.51\pm 0.01$, and the peaks of the γ-ray profiles are at ${\phi }_{1}\,=$ 0.021 ±  0.004 and ${\phi }_{2}\,=$ 0.553 ±  0.013. Adopting this γ-ray peak definition leads to a radio-to-γ-ray lag of $\delta ={\phi }_{1}-{\phi }_{{\rm{r}}}\,\simeq$ 0.04 and a γ-ray peak separation of ${\rm{\Delta }}={\phi }_{2}-{\phi }_{1}\,\simeq$ 0.53. These numbers seem consistent with other LAT MSPs (Abdo et al. 2013); the ${\rm{\Delta }}\gtrsim 0.5$ in phase, however, might indicate that the definition of the first and second γ-ray peaks could, in principle, be reversed. If that is the case, the γ-rays either lead the radio by ∼0.49 or trail it by ∼0.51 in phase.

The alignment of the main peaks of the radio and γ-ray profile might be real, but it could also reflect the limited baseline of the radio timing of the pulsar. A 10⁻³ pc cm⁻³ variation in DM over the length of the Fermi mission could lead to a shift of  ∼0.1 in rotational phase between the radio and the γ-ray profiles. Such a DM variation would be consistent with those seen for other MSPs (Keith et al. 2013). A higher-frequency radio profile (e.g., measured at 820 MHz) would be less sensitive to DM variations, but we have so far been unable to detect PSR J1552+5437 at higher radio frequencies.