LOFAR Discovery of a 23.5 s Radio Pulsar

LOTAAS is an all-northern-sky survey for pulsars and fast transients using LOFAR at a central observing frequency of 135 MHz, with a bandwidth of 32 MHz. The survey employs the six high-band antenna (HBA) stations of the inner core of LOFAR with a maximum baseline of 320 m, known as the Superterp. Each 1 hr pointing consists of three subarray pointings (SAPs), centered at three nearby positions separated by 3.8 degrees. An incoherent beam is formed for each SAP with a full width at half maximum (FWHM) of 55. At the center of each SAP, a hexagonal grid of 61 coherently added tied-array beams (TABs) is formed, each with an FWHM of 21' at the central frequency and overlapping near the half-power point (see Stappers et al. 2011, for more details on LOFAR beamformed modes). A further 12 TABs are placed within each SAP, either toward known sources outside the hexagonal grid, or at other predetermined positions. The Stokes I data are recorded with a sampling time of 492 $\mu {\rm{s}}$ and 2596 channels of 12.2 kHz each.

The raw data are stored in the LOFAR Long Term Archive and then transferred to the Dutch National Supercomputer Cartesius where a presto-based (Ransom 2001; Ransom et al. 2002) pulsar search pipeline is applied. Radio frequency interference (RFI) mitigation is performed on the data using rfifind and then dedispersed using prepsubband to dispersion measures (DMs) from 0 to 550 pc cm⁻³, with step sizes of 0.01 pc cm⁻³ for DMs below 40 pc cm⁻³, 0.05 pc cm⁻³ for DMs between 40 and 130 pc cm⁻³, and 0.1 pc cm⁻³ otherwise. Each dedispersed time series is searched for pulsars with a Fast Fourier Transform (FFT)-based periodicity search with accelsearch, but without searching for accelerated signals. A sifting algorithm that checks whether a candidate is detected in nearby DM values, as well as having a spectral significance higher than 5σ is applied. The sifted candidates are then folded with prepfold, yielding an average of 20,000 folded candidates per pointing. A machine learning classifier (Tan et al. 2018) is then used to determine which of the folded candidates are more likely to be astrophysical in nature. A separate machine learning classifier is also used to search the data for dispersed single pulses (Michilli et al. 2018).

PSR J0250+5854 was discovered in a LOTAAS observation obtained on 2017 July 30. For this observation, the machine learning classifier reported several pulsar candidates at a DM near 45 pc cm⁻³ from a single TAB, with profile significance ranging from 8.4σ down to 4.7σ. The significance of the weakest candidate is similar to some of the most significant candidates automatically classified as nonpulsar. Their spin periods were harmonically related, indicating a fundamental spin period of 23.535 s. We verified that this was the true period by folding the data at this period, shown in Figure 1, as well as several harmonics. Single pulses from the pulsar are also detected in the TAB and these were used to confirm that this was the true period. We formed the time averaged profile of both the odd and even numbered single pulses using this period and the pulsar was seen in both. Folding at half the 23.535 s period, however, resulted in the pulsed signal appearing in only the odd or even numbered average profile.

Zoom In Zoom Out Reset image size

An initial follow-up observation was conducted using the discovery central frequency and bandwidth, but now using all 24 HBA stations of the LOFAR core, providing a longer maximum baseline of 2 km, in order to localize the pulsar. It used a hexagonal grid of 127 TABs centered at the discovery position, each with an FWHM of 38. The pulsar was detected in two adjacent TABs resulting in a signal-to-noise ratio (S/N) weighted position of ${\alpha }_{{\rm{J}}2000}={02}^{{\rm{h}}}{50}^{{\rm{m}}}{11}^{{\rm{s}}}\pm {4}^{{\rm{s}}}$, ${\delta }_{{\rm{J}}2000}=+58^\circ 54^{\prime} \pm 1^{\prime}$. Another pair of follow-up observations centered on this new position, but using a separation between TABs of just 05 resulted in our best initial localization position of ${\alpha }_{{\rm{J}}2000}\,={02}^{{\rm{h}}}{50}^{{\rm{m}}}{17}^{{\rm{s}}}\pm {1}^{{\rm{s}}}$, ${\delta }_{{\rm{J}}2000}=+58^\circ 53^{\prime} 26^{\prime\prime} \pm 11^{\prime\prime}$. However, due to ionospheric beam jitter, these positions have an additional systematic uncertainty of up to 1'.

Inspection of earlier LOTAAS pointings revealed pulsations from PSR J0250+5854 in an observation obtained on 2015 August 5, with the TAB center located 8' from the position of the pulsar. The detection is much weaker than in the discovery observation, as this TAB is located further away from the SAP center compared to the discovery TAB from its SAP center, and as a result has a lower sensitivity. Hence, the pulsar was not detected by the search pipeline in this earlier observation. This fortuitous observation, taken almost exactly two years before the observation with which the pulsar was discovered, minimized any influence of position error on the period determination (see Lorimer & Kramer 2005). This, and the time span between detections, allowed us to initially estimate the spindown rate of the pulsar to be $(3.5\pm 1.4)\times {10}^{-14}$ s s⁻¹.

Ten weekly timing observations of PSR J0250+5854, for an hour per epoch (corresponding to roughly 150 pulse periods) were obtained between 2017 September 9 and November 14 with the HBA antennas of the LOFAR core. Dual polarization, Nyquist sampled time series in complex voltage (CV) mode were recorded for 400 subbands of 195.3 kHz each, centered at 149 MHz. The CV data were combined with the localization data obtained from the LOTAAS pointings and follow-up observations to obtain a timing solution for the pulsar.

A LOFAR Low Band Antenna (LBA) observation was made on 2017 September 23, using all 24 core LBA stations with 300 subbands of 195.3 kHz each, centered at 62.4 MHz. The data were recorded in CV mode.

The LOFAR Two-meter Sky Survey (LoTSS; Shimwell et al. 2017) observed the region of PSR J0250+5854 on 2016 August 31. The 8 hr interferometric imaging observation was processed using a GRID implementation (Mechev et al. 2017) of the standard direction-independent calibration pipeline¹⁷ (van Weeren et al. 2016; Williams et al. 2016). The direction-dependent calibration and imaging was done using the LoTSS-DR1 processing pipeline (T. W. Shimwell et al. 2018, in preparation), which uses KillMS (Tasse 2014; Smirnov & Tasse 2015) to calculate direction-dependent calibration solutions and DDFacet to apply these during the imaging (Tasse et al. 2018). The LoTSS image (Figure 2) shows a single point-like radio source coincident with the best gridded position of PSR J0250+5854.

Zoom In Zoom Out Reset image size

The LoTSS visibilities have sufficient time (1 s subintegrations) and frequency resolution (12.2 kHz channels) to allow gated imaging on the PSR J0250+5854 ephemeris. The ephemeris connecting the discovery, prediscovery, and confirmation observations was used to predict the pulse phase and spin period of PSR J0250+5854 at the time of the LoTSS observation. Correcting for the 10 s dispersion delay across the 120–168 MHz LoTSS band, an on image was made by selecting three 1 s subintegrations around the predicted arrival time of each pulse in the 8 hr integration. The off-image used the remaining subintegrations.

The on–off images (Figure 3) confirm that the point-like radio source is PSR J0250+5854, as the source is absent in the off-image with a 5σ flux upper limit of 0.75 mJy beam⁻¹. In the on-image, PSR J0250+5854 has a total gated flux density of 36.4 mJy. As the on-image uses 12.7% of the 1 s subintegrations of the 8 hr integration, the pulse averaged flux density is 4.6 mJy. The position of PSR J0250+5854 in the on-image is ${\alpha }_{{\rm{J}}2000}={02}^{{\rm{h}}}{50}^{{\rm{m}}}17\buildrel{\rm{s}}\over{.} 781$, ${\delta }_{{\rm{J}}2000}=+58^\circ 54^{\prime} 01\buildrel{\prime\prime}\over{.} 34$, which is consistent with the gridded TAB position described earlier. In agreement with other fields processed using the LoTSS-DR1 pipeline (T. W. Shimwell et al. 2018, in preparation), the astrometric error is 02 and the flux scale is accurate to approximately 20%.

Zoom In Zoom Out Reset image size

The weekly CV data from LOFAR HBA were coherently dedispersed and folded with dspsr (van Straten & Bailes 2011), while the LOTAAS data were incoherently dedispersed and folded. All the folded data were analyzed with psrchive (Hotan et al. 2004). The integrated pulse profiles were referenced against an analytic, noise-free template, obtained by fitting a single von Mises component to the profile with paas, to measure pulse times of arrival (TOAs). A phase-connected timing solution was determined by modeling the TOAs with tempo2 (Edwards et al. 2006; Hobbs et al. 2006). An offset between the TOAs from the LOTAAS and CV data of 623 ± 5 $\mu {\rm{s}}$ is found and corrected by fitting a jump between the two sets of data.¹⁸

We observed PSR J0250+5854 with the GBT on 2017 October 25. The observation was 94 minutes long with 100 MHz of bandwidth at a central observing frequency of 350 MHz. The Stokes I data were recorded with the GUPPI backend (DuPlain et al. 2008).

We observed PSR J0250+5854 with the Lovell Telescope at Jodrell Bank, UK on two separate days; 2017 August 19 and November 30. The 1 hr long observations were recorded using the ROACH backend (Karuppusamy 2011; Bassa et al. 2016) with 384 MHz of bandwidth at a central observing frequency of 1532 MHz.

We observed PSR J0250+5854 with the Nançay Telescope on two separate days; 2018 July 7 and 14 for 43 and 47 minutes, respectively. The observations were recorded using the NUPPI backend (Cognard et al. 2013) with 512 MHz of bandwidth at a central observing frequency of 1484 MHz.

The field of PSR J0250+5854 was observed for ∼600 s with the Position-Sensitive Proportional Counter (PSPC) on board ROSAT during the All-Sky Survey (sequence rs930604n00, data collected between 1990 July and 1991 August; Boller et al. 2016). No X-ray source was detected at the position of the pulsar. Following the prescriptions of Belloni et al. (1994) and using the X-ray image analysis package ximage,¹⁹ we set a 3σ upper limit on the count rate in the 0.1–2.4 keV band of $4.55\times {10}^{-2}$ counts s⁻¹.

Further observations of PSR J0250+5854 were carried out with the Neil Gehrels Swift Observatory's X-Ray Telescope (XRT; Burrows et al. 2005) in Photon Counting mode on 2018 March 23–24 (obs.ID 00010629001, net exposure: 5.8 ks) and 2018 March 28 (obs.ID 00010629002, 4.0 ks). We derived a 3σ upper limit of $2.34\times {10}^{-3}$ counts s⁻¹ in the 0.3–10 keV energy range from the nondetection of the combined data set.

We also inspected the stacked near-uv images of PSR J0250+5854 obtained with the Swift Ultraviolet/Optical Telescope (UVOT). We derived upper limits of 20.25 mag in the U filter (2.2 ks exposure) and 21.21 mag in the UVW2 filter (7.4 ks). As the UVOT data are less constraining than the XRT data, there is no further discussions of these observations below.

Using the LoTSS imaging position, we derived phase-connected timing solutions for PSR J0250+5854 both with and without the inclusion of the 2015 pre-discovery data from LOTAAS. The solutions agreed with each other, indicating that the solution containing the 2015 point accounted for the correct integer number of pulsar rotations since that time. This solution is presented in Table 1, with the position of the pulsar in the P–$\dot{P}$-diagram shown in Figure 4.

Zoom In Zoom Out Reset image size

Table 1.  The Timing Parameters of PSR J0250+5854 Obtained from the Timing Solution Including the Prediscovery TOA from 2015 August 5

Timing Parameters	Values
R.A., ${\alpha }_{{\rm{J}}2000}$	${02}^{{\rm{h}}}{50}^{{\rm{m}}}17\buildrel{\rm{s}}\over{.} 78(3)$
Decl., ${\delta }_{{\rm{J}}2000}$	$58^\circ 54^{\prime} 01\buildrel{\prime\prime}\over{.} 3(2)$
Spin period (s)	23.535378476(1)
Spin period derivative (s s−1)	$2.716(7)\times {10}^{-14}$
Dispersion measure, DM (pc cm−3)	45.281(3)
Epoch of timing solution (MJD)	57973
Solar system ephemeris model	DE405
Clock correction procedure	TT(TAI)
Time units	TCB
Timing span (MJD)	57238.2-58071.9
Number of TOAs	16
Weighted post-fit residual ($\mu {\rm{s}}$)	493
Reduced ${\chi }^{2}$ value	2.9
Derived Parameters
Galactic longitude (degree)	137.8
Galactic latitude (degree)	−0.5
DM distance (kpc)	1.6
Characteristic age (Myr)	13.7
Surface dipole magnetic field strength (G)	$2.6\times {10}^{13}$
Spin-down luminosity (erg s−1)	$8.2\times {10}^{28}$

Note. The parentheses indicate the 1σ uncertainty in the values. The position of the pulsar is fixed to the position obtained from LoTSS. The large reduced ${\chi }^{2}$ value obtained is likely due to each TOA being formed from a limited number of pulses, which could certainly add some jitter compared to the formal uncertainty.

Download table as:  ASCIITypeset image

We estimated the DM of PSR J0250+5854 by measuring the TOA of pulses in five frequency subbands, with central frequencies of 117.7, 133.3, 148.9, 164.6, and 180.2 MHz. The bandwidth of each subband is chosen to be 15.6 MHz, in order to preserve sufficient S/N. To take into account possible profile evolution and scattering across the HBA band (Hassall et al. 2012; Bilous et al. 2016), we first produced five different templates by fitting a single von Mises component to the profile of each subband of a single observation with paas. We then aligned the templates at the point corresponding to half of the peak height on the leading edge and used them to measure the TOA in each subband. The DM of the pulsar is then measured by minimizing the difference in these TOAs with tempo2. The measured DM of 45.281 ± 0.003 pc cm⁻³ gives a distance of approximately 1.6 kpc, using Galactic electron density models from Cordes & Lazio (2002) and Yao et al. (2017).

As there is some evidence for profile evolution we also attempted to refine the DM measurement by modeling the profiles, rededispersed at a DM of 45.281 pc cm⁻³, with two von Mises components. The resulting templates were aligned using two different methods, at the point corresponding to half of the peak height on the leading edge and at the peak of each template. The measured DM values were 45.262 ± 0.003 and 45.304 ± 0.005 pc cm⁻³ respectively. However, visual inspection of the pulse phase versus frequency plot for the lower DM shows some pulse broadening at the lowest frequencies, while at the higher value no significant difference was discernable. We therefore decided to use the data folded at a DM of 45.281 pc cm⁻³ for further analysis.

In order to study the spectrum of PSR J0250+5854, the LOFAR timing observations were flux-calibrated using the method detailed in Kondratiev et al. (2016). Eight flux-calibrated observations were split into four frequency subbands with central frequencies of approximately 119.6, 139.1, 158.7, and 178.2 MHz. The uncertainty in the average flux density of a single frequency subband is conservatively estimated to be 50% (Bilous et al. 2016). The pulsar is also detected by the GBT, allowing us to estimate the flux density at 350 MHz using the radiometer equation (Lorimer & Kramer 2005), with a receiver temperature, T_rec of 23 K, gain of 2 K Jy⁻¹ and effective bandwidth of 60 MHz due to band edge effects and the presence of RFI (Stovall et al. 2014). The sky temperature T_sky in the direction of PSR J0250+5854 is estimated to be 89 K by extrapolating the value of T_sky at 408 MHz (Haslam et al. 1982) with a spectral index of −2.55 (Lawson et al. 1987; Reich & Reich 1988). The uncertainty of this flux density measurement is estimated at 20%. The spectral index of the pulsar is modeled with a power law where ${S}_{\nu }\propto {\nu }^{\alpha }$. The result is shown in Table 2 and Figure 5, with a fitted spectral index of −2.6 ± 0.5 and ${\chi }_{\mathrm{red}}^{2}$ of 0.76.

Zoom In Zoom Out Reset image size

The fitted spectral index suggests that PSR J0250+5854 has a steep spectrum compared to the average pulsar population (e.g., Bates et al. 2013; Bilous et al. 2016). It also agrees with the flux density measurement from LoTSS. The pulsar is not detected with the LOFAR LBA nor the Lovell telescope and Nançay Radio Telescope. We estimated the upper limit on the flux densities using the radiometer equation and the derived system equivalent flux density of 27 kJy in the most sensitive band of the LOFAR LBA between 50 and 60 MHz (van Haarlem et al. 2013) and assuming a detection threshold S/N of 10 and estimated duty cycle of 0.4% based on the measured width of the profile at 129 MHz. The upper limit on the flux density at 1532 MHz from the Lovell telescope is also estimated using the radiometer equation with T_rec of 25 K, T_sky of 5 K, gain of 1 K Jy⁻¹, bandwidth of 384 MHz, estimated duty cycle of 0.3% based on the measured width of the profile at 350 MHz instead, as the pulse width is found to decrease at higher frequencies, and detection threshold S/N of 10. The upper limit on the flux density at 1484 MHz from the Nançay Radio Telescope has T_rec of 35 K, gain of 1.4 K Jy⁻¹ and bandwidth of 512 MHz, while the other values are the same as the 1532 MHz limit. These limits are also shown in Table 2. The upper limits are in slight tension with the spectrum inferred from detections between 100 and 400 MHz. However, RFI strongly affects the LOFAR LBA data, and the high-frequency observations are possibly also affected by scintillation.

Figure 6 shows the integrated pulse profiles of PSR J0250+5854 at central observing frequencies of 350, 168, and 129 MHz, respectively, after correcting for the dispersive delay. The pulsar exhibits significant profile evolution between 350 MHz and the two LOFAR observing frequencies. The profile at 350 MHz shows a double-peaked structure, with a stronger first peak. The pulsar shows similar pulse profiles in the two LOFAR observing bands: a single-peaked structure with a small trailing edge that resembles a scattering tail. However, the similarity between the pulse profiles suggests that this trailing edge is intrinsic to the pulsar instead of due to scattering, unless the frequency-dependent scattering timescale, τ, of the pulsar has a power-law index much smaller than the theoretical predictions of −4 or −4.4 (this has been seen in several pulsars with LOFAR, see Geyer et al. 2017). We measured the profile widths at 129 and 168 MHz by fitting templates of two von Mises functions to account for the trailing edge. The estimated pulse width at half maximum, w₅₀, of the profiles at 129 and 168 MHz, expressed in a fraction of a full rotation, is 0.0035 ± 0.0002 and 0.0034 ± 0.0002 respectively, while the estimated pulse width at 10% of maximum, w₁₀ of both profiles is 0.0075 ± 0.0004. The w₅₀ and w₁₀ of the profile at 350 MHz is estimated to be 0.0026 ± 0.0005 and 0.008 ± 0.001 by fitting the profile with five von Mises functions, with a larger error due to uncertainty in the baseline of the profile.

Zoom In Zoom Out Reset image size

The single-pulse properties of PSR J0250+5854 at LOFAR and GBT observing frequencies show significant differences (see Figure 7). Over eight 1 hr LOFAR observations, the pulsar is shown to switch off on timescales from one rotation up to about 5 minutes (12 rotations). We estimated the nulling fraction of each observation by measuring the pulse energy distribution using penergy and pdist of the psrsalsa suite (Weltevrede 2016). They show two clearly separable Gaussian distributions, corresponding to the pulsar being on and off. We estimated the nulling fraction to range from 9% up to 42% over the separate observations, with an average of 27%. We calculated the modulation index (Weltevrede 2007) of each phase bin of the on-pulse region of a single LOFAR observation, to look for any variation in the single pulses. The modulation indices are found to be consistently around 0.4 across the on-pulse region, which could be due to the nulling. To confirm this, we removed the nulls from the pulse stack and recalculated the modulation indices. We found that only a few phase bins across the profile showed significant modulation, suggesting that there is no other variation in the single pulses apart from the nulling.

Zoom In Zoom Out Reset image size

At 350 MHz there appears to be a steady weak component with occasional highly variable strong pulses seen only in the first component of the integrated pulse profile. To confirm this, we calculated the modulation indices and found that there is indeed strong modulation only in the first component. We then measured the pulse energy distribution of the first component of the single pulses and found a tail of high-energy pulses. We separate 22 of these strong pulses out of the 242 pulses based on the reported energy and produced two integrated profiles, consisting of the strong pulses, and the remainder, dubbed the weak pulses, as shown in Figure 8. These profiles suggest that PSR J0250+5854 has a regular emission where the two components of the profile have roughly equal intensity. However, when the strong pulses are present, the first component becomes more prominent. We are unable to identify any possible nulling at this frequency due to the weak pulses having relatively low S/N.

Zoom In Zoom Out Reset image size

The similarity in rotational parameters between PSR J0250+5854 and XDINSs prompted a search for a possible high-energy counterpart. As the inferred spindown luminosity, $\dot{E}$, is too low for nonthermal emissions, we only considered possible thermal emissions for the X-ray analysis.

The upper limit on the count rate of PSR J0250+5854 from the PSPC is converted to a flux limit by using WebPIMMS,²⁰ assuming a blackbody spectrum with temperature ${kT}=85\,\mathrm{eV}$ (approximately the temperature of the blackbody spectrum of XDINS RX J0720.4–3125; e.g., Viganò et al. 2013) and a DM-derived absorbing column ${N}_{{\rm{H}}}=1.35\times {10}^{21}$ cm⁻² by assuming a 10% ionization fraction in the interstellar medium (He et al. 2013). In the 0.1–2.4 keV band, we obtained absorbed and unabsorbed limiting fluxes of $7.8\times {10}^{-13}$ and $6.5\,\times {10}^{-12}$ erg cm⁻² s⁻¹, corresponding to an upper limit on the luminosity of 2 × 10³³ erg s⁻¹ at 1.6 kpc. This limit is about an order of magnitude higher than the luminosity of RX J0720.4–3125, the brightest XDINS, and therefore does not constitute a severe limit with regard to the hypothesis of an enhanced thermal emission in PSR J0250+5854.

For the nondetection of the dedicated Swift/XRT observation, we explored a wider space of parameters. For the possible blackbody temperatures, we considered the values obtained by fitting the thermal components observed in RX J0720.4–3125 (85 eV), as well as RX J2143.0+0654 (110 eV) and RX J0420.0–5022 (50 eV), the hottest and coolest XDINSs, respectively, and in the "low-magnetic-field" magnetar SGR 0418+5729 (320 eV; the values are from Viganò et al. 2013). For the interstellar absorption, we take a range of values derived from the DM, ${N}_{{\rm{H}}}=1.35\times {10}^{21}$ cm⁻², to that predicted from the H i maps of Dickey & Lockman (1990), ∼9 × 10²¹ cm⁻²; other values derived from extinction and reddening (Schlegel et al. 1998; Schlafly & Finkbeiner 2011) are bounded by these limits. The different limits in the bolometric thermal luminosities, L_bol are shown in Figure 9, together with the expected L_bol from a radiating surface corresponding to a typical neutron star radius of 12 km (Lattimer 2017) at various temperatures.

Zoom In Zoom Out Reset image size

The rotational parameters of PSR J0250+5854 place it at the right end of the P–$\dot{P}$-diagram, shown in Figure 4, where few pulsars are known. Various death lines are presented in the diagram. As the exact mechanism of radio emission is still not well understood, these death lines are based on conditions in the pulsar magnetosphere required for pair cascade production, thought to be essential for the generation of radio emission.

One prevalent model for pair production is the vacuum gap model (Ruderman & Sutherland 1975), where there exists a vacuum gap on the polar cap of a neutron star in which the electric field and magnetic field are nonorthogonal. In order to sustain pair production, the potential difference across this gap must be sufficiently large. As a pulsar spin period increases, the thickness of the gap increases to maintain the potential difference. However, at large spin periods, the increase in thickness is no longer enough to maintain the required potential difference. As a result, pair production and thus radio emission ceases. This period limit depends on the surface dipole magnetic field strength B, which can be inferred from the period derivative, $\dot{P}$, and the curvature radius of the magnetic field lines, r_c. For a purely dipolar field, ${r}_{c}\sim {10}^{8}$ cm would require a potential difference that is too high to produce pairs even in most normal pulsars (see Equation (6) of Chen & Ruderman 1993 for the death line of pulsars with pure dipole magnetic fields). Hence, strong nondipolar fields are generally required near the surface.

Chen & Ruderman (1993) defined a pulsar death valley bounded between the death line for a neutron star with a pure dipole magnetic field (their Equation (6)), and the death line modeled by their Equation (9) (red dashed line in Figure 4) assuming ${r}_{c}\sim R={10}^{6}$ cm, comparable to the radius R of the neutron star, and a polar cap much smaller than possible with a pure dipole field. PSR J0250+5854 is located beyond this death valley, suggesting that an alternative model is needed to explain the emission.

Chen & Ruderman (1993) also proposed a death line for an extreme case of twisted magnetic field lines at the polar cap (their Equation (10)). PSR J0250+5854 has not yet crossed this alternative death line. The presence of such a tangled field has been argued for in the bidrifting behavior of PSR J0815+0939 (Szary & van Leeuwen 2017) by modeling the magnetosphere using a curvature radius ${r}_{c}\sim R={10}^{5}$ cm.

Zhang et al. (2000) addressed the potential drop problem by taking into account relativistic frame-dragging effects (Muslimov & Tsygan 1992; Muslimov & Harding 1997) with a multipole magnetic field configuration. This places the death line (green dashed line, Figure 4) slightly to the left of the location of PSR J0250+5854. Gil & Mitra (2001) argued that the death line for curvature radiation could shift further downward by considering very curved magnetic field lines, with radii of curvature much smaller than the radius of a typical neutron star. Spitkovsky (2006) found a ${(1+{\sin }^{2}\alpha )}^{-1/2}$ relationship between B and the inclination angle α between the rotation and the magnetic axes of a pulsar, which reduces the maximum potential difference across the vacuum gap produced by pulsars of a given spin period by the same factor. This implies that emission from pulsars with large inclination angle would cease at smaller periods compared to those with small inclination angles. Zhang et al. (2000) also considered the case where the radio emission originates from resonant inverse Compton scattering (Zhang et al. 1997; Harding & Muslimov 1998), which places the pulsar death line (green dashed–dotted line) below the location of PSR J0250+5854.

Zhang et al. (2000) also discussed the location of pulsar death lines within the framework of the space-charge-limited flow (SCLF) model (Arons & Scharlemann 1979). The death lines from curvature radiation and inverse Compton scattering (blue dashed and dashed–dotted lines, Figure 4) for this model, in the case of a multipole magnetic field, both lie below the location of PSR J0250+5854. Harding & Muslimov (2011) suggest that the death line for curvature radiation in the SCLF model can be shifted downward with a distorted magnetic field that produces an offset polar cap. Zhou et al. (2017) discussed the effects of different equations of state of a neutron star on the pulsar death line, in which heavier neutron stars could explain the presence of radio pulsars beyond the standard death line.

An alternative interpretation of the process underlying the cessation of pulsar emission is presented in Szary et al. (2014). There, the radio emission, expressed as a fraction of the pulsar spindown energy, is thought to have a maximum possible efficiency. A population synthesis model based on this assumption creates a "death valley" where in the published realization (Figure 5(b) in Szary et al. 2014), long-period pulsars survive significantly past the traditional death lines. Using 1400 MHz luminosities, a radio efficiency upper limit of 0.01 is derived: For PSR J0250+5854, we have determined a 10σ upper limit on the 1400 MHz period-averaged flux density from the Lovell telescope of 15 μJy. Using the 1.6 kpc distance estimate from DM measurement, the radio luminosity as defined in Equation (2) of Szary et al. (2014) is $3.0\,\times {10}^{26}$ erg s⁻¹. A comparison with the spindown energy of $9\times {10}^{28}$ erg s⁻¹ (Table 1) produces a required efficiency for PSR J0250+5854 of <0.004, below the 0.01 upper bound proposed by Szary et al. (2014) for successful generation of radio emission. This suggests that indeed, dim, long-period pulsars may continue to shine.

Nevertheless, the more recent death line models discussed above are able to explain the presence of PSR J0250+5854. While the pulsar is not as constraining as the 8.5 s PSR J2144−3933 due to its large $\dot{P}$, it still rules out the conventional model described by Equation (9) of Chen & Ruderman (1993). This also suggests that there could be a lot more pulsars in this region of the P–$\dot{P}$-diagram that are yet to be found.

The similarity of the rotational parameters of PSR J0250+5854 to the XDINSs and magnetars indicates a possible connection between them. Viganò et al. (2013) suggest that pulsars with rotational parameters similar to those of PSR J0250+5854 could evolve from magnetars with an initial dipolar magnetic field strength on the order of 10¹⁵G with a large decay in magnetic field. Using their Figure 10, we estimated the real age of PSR J0250+5854 to be ∼3 × 10⁵ years, much smaller than the characteristic age, with an expected thermal luminosity of 10³⁴ erg s⁻¹ extrapolated from their Figure 11. However, the Swift/XRT nondetection suggests that such an evolutionary track is unlikely, as it places the upper limit on ${kT}=85\,\mathrm{eV}$ for ${N}_{{\rm{H}}}=9\times {10}^{21}$ cm⁻², which would require an emitting surface with a radius three times larger than that of a canonical neutron star. PSR J0250+5854 does not show any evidence of a magnetar nature at the moment. Though transient, radio-loud magnetars have very different properties compared to PSR J0250+5854 (e.g., their spindown luminosity is much larger and their characteristic age is much lower), it is still possible, however unlikely, that it possesses a strong toroidal component and will manifest itself as a magnetar in the future.

The connection between PSR J0250+5854 and XDINSs is still possible if XDINSs have a nonmagnetar origin. The nondetection does not place a hard constraint as most of the known XDINSs have a lower L_bol than our most optimistic upper limit of ${L}_{\mathrm{bol}}=8\times {10}^{31}$ erg s⁻¹ for ${kT}=110\,\mathrm{eV}$ and ${N}_{{\rm{H}}}=1.35\times {10}^{21}$ cm⁻². However, it would be difficult to detect soft X-ray emission (0.1–2.4 keV) from PSR J0250+5854 due to the large ${N}_{{\rm{H}}}$ along the line of sight to the pulsar, which absorbs a larger fraction of the soft X-ray flux.

Deep X-ray and optical observations of the 8.5 s radio pulsar PSR J2144−3933 showed no evidence of enhanced thermal emission (L_bol > $\dot{E};$ Tiengo et al. 2011). It is possible that PSR J0250+5854 also shows no such emission based on the nondetections.

The small duty cycle of PSR J0250+5854 is expected, given that the pulsar beam radius (for nonrecycled pulsars), ρ, appears to scale as $\rho =K/\sqrt{P}$ (see, e.g., Lorimer & Kramer 2005 and references therein). The factor K depends weakly on frequency. At 400 MHz, $K\sim 8$ deg ${{\rm{s}}}^{1/2}$ for a beam radius measured at 10% of the peak intensity, ${\rho }_{10}$ (e.g., Mitra & Rankin 2002). For this pulsar, we expect ${\rho }_{10}\sim 1\buildrel{\circ}\over{.} 7$. At a 10% intensity level of the outer pulse edges, we measure a width of ${w}_{10}\sim 3^\circ$. Therefore, within the uncertainties, ${w}_{10}\sim 2{\rho }_{10}$, and hence we can confirm the general expected scaling, which is based on the notion that with increasing light cylinder radius, the open field-line regions narrow, with a predominantly dipolar magnetic field structure as indicated by the $\sqrt{P}$ relationship. Note that the observed pulse width is determined at the height where radio emission leaves the magnetosphere, which is typically on the order of hundreds of kilometers (Mitra & Rankin 2002; Kijak & Gil 2003). This is far above the neutron star surface where pair production takes place. While the multipole components of the field are required to push the pulsar death line to longer periods (as argued in Section 4.1), they are not necessarily significant anymore at the emission height given their rapid decay with distance. Hence at the emission height the magnetosphere is approximately dipolar.

The exact dependence between w₁₀ and ${\rho }_{10}$ depends on the geometrical configuration, i.e., the magnetic inclination angle α and the impact angle β (see, e.g., Lorimer & Kramer 2005). Without polarization data measured over a sufficiently large duty cycle, it is difficult to determine these angles (e.g., from the position angle swing in terms of a rotating vector model; Radhakrishnan & Cooke 1969). Occasionally, it is possible to use statistical arguments to derive constraints (e.g., Kramer et al. 1998), but we have failed in this particular case, as a wide range of $\alpha -\beta$ combinations are consistent with the observed width and inferred beam radius, given, in particular, the uncertainty in the K factor when scaling to such large periods.

Nevertheless, it is clear that overall the beam radius is small. The corresponding beaming fraction, f (i.e., the fraction of celestial sphere covered by the sweeping beam) is also small, reducing the likelihood to detect such long-period pulsars for randomly oriented beam directions. Assuming a random distribution in α between 0° and 180°, one finds an expectation value for $f=(1-\cos \rho )+(\pi /2-\rho )\sin \rho$ (Emmering & Chevalier 1989), where ρ is expressed in radians. Inserting our inferred ${\rho }_{10}$, this implies $f\sim 0.05$, so only 5% of such long-period pulsars are observable, and the number of potentially existing pulsars is correspondingly larger. However, if α is small, for instance, due to an alignment of the pulsar magnetic axis with age (e.g., Johnston & Karastergiou 2017), then the beaming fraction is much smaller still, and even more such pulsars could exist. Other effects may also alter the effective beam radius, such as unfilled beams (Lyne & Manchester 1988) or different emission heights as seen in different pulsars (Lorimer & Kramer 2005). Overall, we find it remarkable that the observed pulse width essentially follows the extrapolation from smaller periods.

The period dependence of f as shown by PSR J0250+5854 is important in pulsar population synthesis for explaining why pulsars do not pile up near the death line (Hartman et al. 1997; van Leeuwen & Stappers 2010; Johnston & Karastergiou 2017), without needing to invoke an alternative explanation of magnetic field decay in nonrecycled, isolated pulsars (see van Leeuwen & Verbunt 2004; Johnston & Karastergiou 2017). Furthermore, the probable nonmagnetar origin (as discussed in Section 4.2) and determination of a $3\times {10}^{13}$ G field in a 23.5 s pulsar (Table 1) argues against magnetic field decay.

Our inability to obtain a definitive measurement of DM for PSR J0250+5854 highlights the difficulty in disentangling the pulse shape evolution at the large fractional bandwidth and low frequencies of LOFAR (see Hassall et al. 2012). The different templates produced for each of the subbands of the pulsar suggest that there is profile evolution across the bandwidth of the LOFAR observation. However, the evolution does not look like scattering as the width of the templates does not increase at lower frequencies and in fact we seem to lose a component seen at higher frequencies. The frequency evolution of the profiles also means that there is no obvious reference point in aligning the templates to measure the DM. For example, when we align the templates at the peaks, the TOAs generated at 180.2 MHz do not line up with others after fitting for DM.

We have shown that more pulsars with similarly long periods as PSR J0250+5854 could exist, and if they are detectable, will likely have small duty cycles as well. This would make FFT-based periodicity searches less sensitive to these sources. Hence, many of the longer period pulsars discovered recently by various surveys (Keane et al. 2010; Karako-Argaman et al. 2015; Deneva et al. 2016) were first detected via single pulses as RRATs, with the underlying periodicity only detected with longer follow-up observations. However, PSR J0250+5854 was discovered through an FFT-based periodicity search. This is primarily due to the relatively consistent flux from the pulsar apart from the nulling, combined with the 1 hr long observations of LOTAAS giving us a large number of pulses (∼150) for detection. For comparison, the two-minute dwell time of the GBNCC survey (Stovall et al. 2014) would only see 5–6 pulses. It should be noted that the pulsar was only detected at the fifth and higher harmonics in the discovery observation. While the fifth harmonic has a detection significance of 9σ as reported by accelsearch, the lower harmonics are not detected down to significance of 2σ. The FFT-based periodicity search of LOTAAS is restricted to a minimum Fourier domain frequency of 1/16 Hz and a maximum number of harmonics of 16. This meant that the fundamental frequency of the pulsar will not be detected by the search pipeline. The lower harmonics of the pulsar were not detected partly due to the presence of red noise in the low-frequency regime. The small duty cycle of the pulsar also meant that the power of the pulsar in the Fourier domain is spread over a large number of harmonics, thus requiring the harmonic summing in the Fourier spectrum of a large number of harmonics for the pulsar to be detectable (Manchester & Taylor 1977). Fortunately, the large flux density of the pulsar at LOTAAS observing frequencies resulted in large power of the higher harmonics, allowing us to detect the pulsar through the summation of higher harmonics.

The discovery of PSR J0250+5854 was also aided by the fact that the S/N shows a clear and narrow peak in DM space (the peak is only 4 pc cm⁻¹ wide). The pulsar is thus well distinguished from zero-DM RFI, which would be more difficult at higher observing frequencies. The discovery was also helped by the steep spectrum of PSR J0250+5854. If other long-period pulsars showed similar spectral indices then they would be easier to detect with LOFAR as well.

Our observation of PSR J0250+5854 at a frequency of 350 MHz showed somewhat different emission properties compared to the LOFAR observations. The pulsar showed sporadic strong single pulses at 350 MHz. This suggests that the emission mechanism is strongly frequency dependent and perhaps results in less stable emission at higher frequencies. It is possible that there are other long-period pulsars that show variable single pulses at higher observing frequencies, which are more likely to be discovered through single pulse searches. On the other hand, they might exhibit more regular emission at lower frequencies that would require periodicity searches and a long dwell time that LOTAAS provided to be detected.

The discovery of PSR J0250+5854 has significantly expanded the known range of rotation-powered pulsar periods. However, FFT-based periodicity searches will likely miss them if they have lower flux compared to PSR J0250+5854. The presence of red noise would also reduce the sensitivity toward long-period pulsars. Furthermore, if PSR J0250+5854 had a larger duty cycle, the pulsar would not have been detectable, due to the smaller number of significant harmonics. Hence it is likely that we are missing out on more long-period pulsars due to these shortcomings of FFT-based periodicity searches. Recently, there has been renewed interest in using the Fast Folding Algorithm (FFA, Staelin 1969) to search for pulsars. Various different implementations of FFA have been developed or are currently in preparation (Kondratiev et al. 2009; Cameron et al. 2017; V. Morello et al. 2018, in preparation; Parent et al. 2018) with the expectation that FFA will be more sensitive toward long-period pulsars. In fact, in a preliminary analysis, the FFA implementation of V. Morello et al. (2018, in preparation) detected PSR J0250+5854 at the fundamental period with a high S/N of 40, much higher than the significance of the detection made by the FFT. The addition of FFA to the processing pipeline of LOTAAS is ongoing, and it will potentially discover new pulsars that will populate the area in the P–$\dot{P}$-diagram around PSR J0250+5854.