LOFAR Discovery of the Fastest-spinning Millisecond Pulsar in the Galactic Field

The discovery of the first millisecond pulsar (MSP), PSR B1937+21 with a spin frequency of 642 Hz, by Backer et al. (1982) came as a great surprise, and demonstrated that some neutron stars can reach astounding rotational rates. The more recent discovery of transitional millisecond pulsars (tMSPs), which transition back and forth between a rotation-powered MSP and an accretion-powered low-mass X-ray binary (LMXB) state (Archibald et al. 2009; Papitto et al. 2013; Bassa et al. 2014), confirmed the basic recycling model of Alpar et al. (1982) and Radhakrishnan & Srinivasan (1982) in which a neutron star is spun up to millisecond spin periods due to the accretion of matter and angular momentum. At the same time, the tMSPs have also raised many questions about the detailed physics of the pulsar recycling process and how efficient it can ultimately be in terms of spinning up neutron stars (e.g., Archibald et al. 2015; Deller et al. 2015; Papitto et al. 2015; Jaodand et al. 2016).

It is striking that, since the discovery of PSR B1937+21, only one faster-spinning MSP has been found (PSR J1748−2446ad, in the globular cluster Terzan 5, spinning at 716 Hz; Hessels et al. 2006). While the neutron star equation-of-state in principle allows spin frequencies up to 1200 Hz (Cook et al. 1994; Lattimer & Prakash 2004) before mass-shedding or break-up, the currently observed spin frequency distribution of radio and X-ray MSPs cuts off around 730 Hz (Chakrabarty et al. 2003; Ferrario & Wickramasinghe 2007; Chakrabarty 2008; Hessels 2008). The question thus remains: can nature form submillisecond pulsars?

Physical effects such as decoupling of the Roche lobe (Tauris 2012), transient accretion (Bhattacharyya & Chakrabarty 2017), and gravitational-wave emission (Chakrabarty et al. 2003) have been put forward to explain the observed spin frequency cutoff. Observationally, there are additional challenges in detecting submillisecond pulsars, compared to canonical MSPs with spin frequencies between 200 and 500 Hz, but the cutoff around 730 Hz is hard to explain purely as an observational bias. For example, computational advances allow present-day radio pulsation surveys to retain sensitivity to spin frequencies well in excess of 1000 Hz (e.g., Lazarus et al. 2015) because it is now possible to record data with sufficient time and frequency resolution—which is critical for correcting for the dispersive delays introduced by the ionized interstellar medium (IISM). Furthermore, if MSP searches are conducted at sufficiently high radio frequencies (1–2 GHz) the effects of scattering in IISM should also not preclude the detection of submillisecond pulsars, though high-frequency radio searches are disadvantaged by the fact that pulsars typically have steep radio spectra ($S\propto {\nu }^{\alpha }$, where $\alpha =-1.4\pm 1.0;$ Bates et al. 2013).

Watts et al. (2015) suggest that a possible bias against finding rapidly spinning pulsars, and hence energetic MSPs, may be the irradiation-driven mass loss from the binary companion, which can lead to eclipses of the radio signal during large parts of the orbit (Stappers et al. 2014; Roy et al. 2015). Similar considerations were previously presented by Tavani (1991), motivated by the discovery of the first eclipsing black widow pulsar binary system, PSR B1957+20, in which the low-mass, bloated companion star is irradiated by the pulsar wind (Fruchter et al. 1988). Indeed, there is evidence that the eclipsing MSP systems—both the black widows with very-low-mass companions and the redbacks with higher-mass (${M}_{{\rm{c}}}\gtrsim 0.2$ ${M}_{\odot }$) companions—are on average spinning faster than "classical" MSPs with white dwarf companions (Hessels 2008; Papitto et al. 2014). For eclipsing systems, there is again an advantage toward observing at higher radio frequencies, where the eclipse durations are typically lower (Archibald et al. 2009), but also the disadvantage that the intrinsic pulsar spectrum is generally falling off rapidly toward higher frequencies.

Recent results by Kuniyoshi et al. (2015), Kondratiev et al. (2016), and Frail et al. (2016) indicate that the fastest-spinning MSPs tend to have the steepest radio spectra ($\alpha \lt -2.5$), pointing to another possible bias against finding fast-spinning MSPs in ongoing surveys, which focus on central observing frequencies around 350 MHz (Deneva et al. 2013; Stovall et al. 2014) and 1.4 GHz (Cordes et al. 2006; Keith et al. 2010; Barr et al. 2013). As a result, radio pulsation searches at frequencies below 300 MHz have the potential of opening up a so far largely unexplored parameter space, in the cases where IISM scattering is low and eclipsing does not hinder detection either.

Here, we present the discovery of PSR J0952−0607, a very-steep-spectrum MSP, which is now the fastest-spinning neutron star known in the Galactic field (outside of a globular cluster). PSR J0952−0607 was found in a radio pulsation survey using the Low-Frequency Array (LOFAR; Stappers et al. 2011; van Haarlem et al. 2013) to target unassociated Fermi γ-ray sources. This Fermi-targeted approach has been successful in finding many new MSPs (Ray et al. 2012), but our survey is the first to use LOFAR to survey at observing frequencies of 135 MHz (see also Pleunis et al. 2017). To enable this survey, a combination of coherent and incoherent dedispersion has been employed to limit the effects of dispersive smearing (Bassa et al. 2017). In Section 2, we will highlight the discovery and multi-wavelength follow-up, with the results being presented in Section 3. We discuss the broader implications of PSR J0952−0607's discovery and conclude this manuscript in Section 4.

2.1. Radio

PSR J0952−0607 was discovered as part of an ongoing LOFAR survey at 135 MHz, continuing on the pilot survey by Pleunis et al. (2017). Unassociated γ-ray sources were selected from an all-sky source list based on seven years of Fermi-Large Area Telescope (LAT; Atwood et al. 2009) Pass 8 data. That list resulted from a preliminary version of the procedure that will be used to produce the next public release LAT source catalog. Among these is a new γ-ray source that has a test statistic of $\,100$ and a relatively small error radius of $3\buildrel{\,\prime}\over{.} 8$. The source's γ-ray spectrum is strongly curved, peaking at 1.4 GeV, and no significant γ-ray emission is detected above 10 GeV. As the source lies well out of the plane ($b=35\buildrel{\circ}\over{.} 4$), these characteristics made it a prime millisecond pulsar candidate.

The γ-ray source was observed for 20 minutes on 2016 December 25 with LOFAR. The high-band antennas (HBAs) from the innermost 21 LOFAR core stations (longest baseline of 2.3 km) were used to form seven tied-array beams ($3\buildrel{\,\prime}\over{.} 5$ FWHM), each covering a 39 MHz wide band centered at 135 MHz. This setup is identical to that of Pleunis et al. (2017). The complex voltage output for each tied-array beam—in the form of 200 Nyquist-sampled, dual-polarization subbands of 195 kHz each—was processed using GPU-accelerated software to perform coherent dedispersion and channelization with cdmt (Bassa et al. 2017) at steps of 1 pc cm⁻³ between dispersion measures (DMs) of 0.5 and 79.5 pc cm⁻³. The resulting coherent filterbanks, sampled at 81.92 μs and 48.83 kHz in time and frequency, were dedispersed incoherently around the coherently dedispersed DM trial (${\rm{\Delta }}\mathrm{DM}=-0.5$ to 0.5 pc cm⁻³) at steps of 0.002 pc cm⁻³ using the Dedisp library (Barsdell et al. 2012). The dedispersed time series were searched for periodic signals using frequency domain acceleration searching (tools from PRESTO; Ransom 2001; Ransom et al. 2002). A detailed description of the processing steps is given in Bassa et al. (2017).

PSR J0952−0607 was discovered blindly at high significance in four of the seven tied-array beams at a spin frequency of 707 Hz and a DM of 22.41 pc cm⁻³. The pulsar was found at an acceleration of 1.3 m s⁻², indicating that it is part of a binary system. As the cumulative pulse profile of PSR J0952−0607 is double peaked, with components separated by approximately 110° (see Figure 1), we verified that the 707 Hz spin frequency is the fundamental (i.e., the neutron star's true rotation rate) by folding the dedispersed time series of the discovery observation at several harmonically related spin frequencies. In all cases, the resulting profiles were the sum of copies of the 707 Hz profile and were of lower signal-to-noise ratio than the 707 Hz profile.

Zoom In Zoom Out Reset image size

Follow-up observations (10 minute integration times) were obtained with the HBAs from 23 LOFAR core stations (longest baseline of 3.5 km) using 7 tied-array beams with 39 MHz of bandwidth centered at 135 MHz on 2017 January 4 (initial follow-up gridding observation) and a single beam with the full HBA band (78 MHz at 149 MHz) for all subsequent observations. These observations allowed us to refine the position of the pulsar and start the timing program. A 3 hr HBA integration was obtained on 2017 January 28/29 to constrain the orbital parameters. To determine the radio spectrum of PSR J0952−0607, we obtained a 2 hr integration with the LOFAR low-band antennas (LBAs) between 30 and 90 MHz on 2017 February 5/6 and a 47 minute observation at 350 MHz (100 MHz bandwidth) on 2017 March 1 with GUPPI (DuPlain et al. 2008) at the Green Bank Telescope (GBT). The pulsar was not detected in the LOFAR LBA observation, but was easily seen in the GBT 350 MHz observation.

The complex voltage data of the discovery and follow-up LOFAR HBA observations were coherently dedispersed and folded with dspsr (van Straten & Bailes 2011) and analyzed using psrchive (Hotan et al. 2004) tools. Pulse profiles for two-minute sub-integrations were referenced against an analytical pulse profile template to obtain time-of-arrival (TOA) measurements. A phase-connected timing solution, accounting for every rotation of the pulsar, was determined from these TOAs using tempo2 (Edwards et al. 2006; Hobbs et al. 2006).

2.2. Optical

2.3. X-Ray

The phase-connected timing solution models the rotation and orbit of PSR J0952−0607. As the timing solution has a time baseline of approximately a third of a year, the spin parameters (ν and $\dot{\nu }$) are degenerate with the astrometric parameters (${\alpha }_{{\rm{J}}2000}$, ${\delta }_{{\rm{J}}2000}$). Assuming values for the spin frequency derivative $\dot{\nu }$ between $-{10}^{-12}$ and $-{10}^{-15}$ s⁻², the fitted position of the pulsar traces a line on the sky as indicated in Figure 2. Along this line the optical images show a strongly variable object, located about $1^{\prime}$ from the LOFAR-gridded position of PSR J0952−0607. The object varied by at least 1.5 mag in the co-added ${r}^{{\prime} }$ images, being below the detection threshold in approximately half of them. The ${r}^{{\prime} }$-band magnitudes are modulated at the orbital period of PSR J0952−0607 (Figure 3), confirming that the object is the binary companion of the pulsar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The binary companion to PSR J0952−0607 is located at ${\alpha }_{{\rm{J}}2000}={09}^{{\rm{h}}}{52}^{{\rm{m}}}08\buildrel{\rm{s}}\over{.} 319(3)$ and ${\delta }_{{\rm{J}}2000}=-06^\circ 07^{\prime} 23\buildrel{\prime\prime}\over{.} 49(5)$. The positional uncertainty quoted is the quadratic sum of the uncertainty in the astrometric calibration and the positional uncertainty of the companion on the co-added image (of order $0\buildrel{\prime\prime}\over{.} 05$). To estimate the impact of the uncertainty in the optical position of the binary companion on the parameters of the timing solution, we performed a Monte Carlo simulation, drawing positions (${\alpha }_{{\rm{J}}2000}$, ${\delta }_{{\rm{J}}2000}$) from normal distributions with appropriate means and widths. For each of these positions, the remaining parameters in the timing solution were fitted to generate distributions from which the parameters and their uncertainties were determined. These values are listed in Table 1. Taking into account the positional uncertainties yields a 3σ limit on the spin frequency derivative $\dot{\nu }\gt -3.3\times {10}^{-15}$ s⁻², corresponding to a spin period derivative of $\dot{P}\lt 1.1\times {10}^{-20}$ s s⁻¹, placing it in the lower half of the MSP $\dot{P}$ distribution. Because of the short spin period, the surface magnetic field $B\propto \sqrt{P\dot{P}}$ is low, $B\lt 1.3\times {10}^{8}$ G. The limit on the spin-down luminosity ($\dot{E}\propto \dot{P}{P}^{-3}$) is $\dot{E}\lt 1.6\times {10}^{35}$ erg s⁻¹. We note that these values could be significantly lower should PSR J0952–0607 have a significant proper motion, as the observed $\dot{P}$ would have a contribution due to the Shklovskii effect (Shklovskii 1970).

Table 1.  Parameters for PSR J0952−0607

Parameters	Value
Timing Parameters
R.A., ${\alpha }_{{\rm{J}}2000}$	${09}^{{\rm{h}}}{52}^{{\rm{m}}}08\buildrel{\rm{s}}\over{.} 319(3)$
Decl., ${\delta }_{{\rm{J}}2000}$	$-06^\circ 07^{\prime} 23\buildrel{\prime\prime}\over{.} 49(5)$
Spin frequency, ν (s−1)	$707.314434911(16)$
Spin frequency derivative, $\dot{\nu }$ (s−2)	$\gt -3.3\times {10}^{-15}$
Epoch of timing solution (MJD)	57800
Dispersion measure, DM (pc cm−3)	$22.41149(10)$
Binary model	ELL1
Orbital period, Pb (day)	$0.267461038(12)$
Projected semimajor axis, x (s)	$0.0626694(14)$
Time of ascending node passage, Tasc (MJD)	$57799.9119800(8)$
Solar system ephemeris model	DE421
Clock correction procedure	TT(BIPM2011)
Time Units	TCB
Timing Span (MJD)	57747.1–57851.9
Number of TOAs	164
Weighted rms post-fit residual (μs)	5.6
Reduced ${\chi }^{2}$ value	1.11

Note. The astrometric parameters (${\alpha }_{{\rm{J}}2000}$ and ${\delta }_{{\rm{J}}2000}$) are kept fixed at the position of the optical counterpart. The eccentricity is kept fixed at e = 0, implicitly assuming that the orbit is circular. For the ELL1 binary model (Lange et al. 2001), this means $\kappa =e\sin \omega =0$ and $\eta =e\cos \omega =0$.

Download table as:  ASCIITypeset image

Models for the electron distribution in our Galaxy along the line-of-sight to PSR J0952−0607 ($l=243\buildrel{\circ}\over{.} 65$, $b=35\buildrel{\circ}\over{.} 38$) constrain the distance through the observed DM. The NE2001 model by Cordes & Lazio (2002) predicts a distance of $d=0.97\,\mathrm{kpc}$, while the YMW16 model (Yao et al. 2017) places it significantly farther away at $d=1.74\,\mathrm{kpc}$. At either distance, the Galactic extinction model by Green et al. (2014, 2015) estimates the same reddening of ${E}_{B-V}=0.061$, leading to an extinction of ${A}_{{r}^{{\prime} }}=0.14$ (Schlafly & Finkbeiner 2011).

The orbital parameters from the timing solution (Table 1) show that PSR J0952−0607 is in a 6.42 hr binary around an ultra-low-mass companion. The mass function is $3.69\times {10}^{-6}$ ${M}_{\odot }$, setting the minimum companion mass at 0.019 M${}_{\odot }$ for a 1.4 M${}_{\odot }$ pulsar. These properties are consistent with PSR J0952−0607 being a black widow system, where matter from the companion is ablated by the energetic pulsar wind (e.g., Fruchter et al. 1988; see Roberts 2013 for a review).

The black widow nature of PSR J0952−0607 is confirmed by the sinusoidal optical light curve, which is consistent with irradiation of the companion hemisphere facing the pulsar. We use the Icarus software (Breton et al. 2012) to model the ${r}^{{\prime} }$-band light curve. As the absence of color information precludes a full parameter fit, we make the following assumptions to estimate system parameters. We assume that the companion is co-rotating and that its base temperature (that of the unirradiated hemisphere) is 2500 K—in line with what is seen in other black widows (see, e.g., Stappers et al. 2001; van Kerkwijk et al. 2011; Romani et al. 2016)—and hence contributes minimally to the flux of the irradiated hemisphere. The reddening is kept at ${E}_{B-V}=0.061$ (see above), and the pulsar mass is fixed at 1.4 M_⊙. Finally, we require the dayside temperature to be such that the implied irradiation represents of order 10% of the spin-down luminosity for an assumed credible range from 5 × 10³⁴ to $1.6\times {10}^{35}$ erg s⁻¹ (Breton et al. 2013).

We find that the observed light curve is mostly inconsistent with models that assume a filling factor of unity (i.e., where the companion fills the Roche lobe), as it leaves the orbital inclination largely unconstrained, though leaning toward edge-on ($i=90^\circ$), with tightly correlated dayside temperatures and distances in the range of 4500–5800 K and 2.3–5.3 kpc, respectively. The goodness of fit improves significantly for models with an assumed filling factor of 0.5, providing a well-defined inclination of $i\sim 40^\circ$. The correlation between dayside temperature and distance is weaker, yielding similar dayside temperatures and slightly smaller distances (1.7–3.8 kpc). Models with the filling factor as a free parameter prefer slightly smaller filling factors while the overall goodness of the fit is not significantly increased. These models provide similar constraints on the inclination, dayside temperature, and distance, with even weaker correlations due to the extra free parameter.

We conclude that PSR J0952−0607 is almost certainly not close to Roche-lobe filling, with an orbital inclination in the intermediate range. The predicted model distances are at the high end of those estimated from the DM, which may suggest Roche-lobe filling factors of 0.5 or less, or lower dayside temperatures, indicating spin-down luminosities of a few 10³⁴ erg s⁻¹ or that the conversion of spin-down into heating is less than 10% efficient. Upcoming multi-color photometry will be able to confirm these values.

The low Roche-lobe filling factor and low inclination are consistent with the absence of eclipses of the radio signal in the observations obtained so far. Eclipses are seen in the majority of black widow systems, and eclipses are generally most pronounced at low observing frequencies (e.g., Stappers et al. 1996; Archibald et al. 2009). Besides eclipses, ionized matter passing through the line-of-sight leads to increases in the DM near orbital phase $\phi =0.25$, resulting in delays in the TOAs. The TOA residuals in Figure 3 do show systematic delays at $\phi =0.26$–0.28, possibly hinting at the presence of ionized material in the line-of-sight. The ongoing LOFAR timing observations will constrain whether material is ablated from the companion of PSR J0952−0607.

We use the Hamaker (2006) model for the LOFAR HBA beam together with the radiometer equation and the method as detailed in Kondratiev et al. (2016) to obtain flux density measurements for PSR J0952−0607 over the HBA band between 110 and 188 MHz. Measurements from five different observations, totaling 1 hr of integration time, were averaged. The flux density at 350 MHz was measured from the single GBT observation, using the radiometer equation with gain, system temperature, and bandwidth values from Stovall et al. (2014). We obtain ${S}_{\mathrm{mean}}=45$, 32, 21, 9, and 1.5 mJy at frequencies of 119.6, 139.2, 158.7, 178.2, and 350 MHz. Following Bilous et al. (2016), we conservatively estimate 50% uncertainties in the flux density measurements. Modeling the spectrum with a power law ${S}_{\nu }\propto {\nu }^{\alpha }$ over frequency ν yields a spectral index of $\alpha =-3.3\pm 0.3$ and a flux density at 150 MHz of ${S}_{150}=21\pm 2$ mJy. We note that PSR J0952−0607 is not detected in the 150 MHz TGSS-ADR source catalog (Intema et al. 2017), listing sources brighter than $7\sigma$ significance with a median noise of 3.5 mJy beam⁻¹. At the location of PSR J0952−0607 the flux density in the TGSS-ADR images is 12 mJy beam⁻¹, suggesting a $3\sigma$ detection, within $2\sigma$ from the LOFAR flux density. Even if the LOFAR fluxes are overestimated by a factor of two, as suggested by Frail et al. (2016), the spectral index remains as steep as $\alpha =-2.6\pm 0.4$.

Figure 1 shows the cumulative pulse profile of PSR J0952−0607 at different frequencies, revealing significant evolution with observing frequency. The pulse profile evolves from two double-peaked components at an observing frequency of 350 MHz to two scatter-broadened components at LOFAR frequencies. To estimate the scattering timescale, we approximate scatter broadening as a convolution with a truncated exponential, appropriate for a thin scattering screen (Williamson 1972), and assume that scatter broadening can be neglected in the 350 MHz pulse profile such that it can be treated as the intrinsic profile. The four components are modeled with von Mises functions and the amplitude of the four components was allowed to vary with frequency, while positions and widths are kept fixed. Using this model we find scattering times of 47 and 113 μs at observing frequencies of 168 and 129 MHz, respectively. The first two components of the pulse profile, measured against the highest component, increase by a factor of 1.9 and 2.0 from 350 to 129 MHz, while the third component decreases by a factor of 3.3.

These scattering timescales indicate that the scatter broadening exceeds the pulse period for observing frequencies below 70 MHz, assuming scattering scales as $\tau \propto {\nu }^{-4}$. As the LBA sensitivity peaks near 60 MHz, where the pulsed signal will be scattered out completely, it is not surprising that PSR J0952−0607 is not detected with the LOFAR LBA.

The Swift/XRT X-ray non-detection of PSR J0952−0607 translates to a $3\sigma$ flux limit in the 0.3–10 keV band of ${f}_{{\rm{X}}}\lt 1.1\times {10}^{-13}$ erg s⁻¹ cm⁻² for absorbed blackbody (${T}_{\mathrm{eff}}=0.23$ keV) and power-law (${\rm{\Gamma }}=2$) spectra. Here, we assumed ${N}_{{\rm{H}}}=4\times {10}^{20}$ cm⁻² estimated from the optical reddening through the relation by Güver & Özel (2009). The resulting $3\sigma$ X-ray luminosity limits (0.3–10 keV) are ${L}_{{\rm{X}}}\lt 1.1\times {10}^{31}$ erg s⁻¹ at a distance of 0.97 kpc and ${L}_{{\rm{X}}}\lt 3.6\times {10}^{31}$ erg s⁻¹ at 1.74 kpc. These limits are consistent with the observed relation between the X-ray luminosity and spin-down luminosity of radio MSPs (Possenti et al. 2002).

PSR J0952−0607 has a spin frequency $\nu =707\,\mathrm{Hz}$. This makes it the fastest-spinning neutron star known in the Galactic field (outside of a globular cluster), surpassing the 35 year record set by the first MSP to be discovered, PSR B1937+21, which spins at 642 Hz (Backer et al. 1982). Only PSR J1748−2446ad, located in the globular cluster Terzan 5, spins faster at 716 Hz (Hessels et al. 2006). Of the 213 known Galactic field MSPs with $P\lt 30$ ms ($\nu \gt 33$ Hz) from the pulsar catalog¹⁹ (Manchester et al. 2005), only 13 have $P\lt 2$ ms ($\nu \gt 500$ Hz). A further 3 MSPs in globular clusters also satisfy this condition, including PSR J1748−2446ad. Of the accreting millisecond X-ray pulsars, only 3 out of 15 have spin frequencies above 500 Hz (Patruno & Watts 2012).

For PSR J1748−2446ad it is not possible to determine the intrinsic $\dot{\nu }$ of the neutron star because the observed change in spin rate with time is dominated by acceleration in the gravitational potential of Terzan 5 (Prager et al. 2016). Conversely, it will be possible to measure PSR J0952−0607's intrinsic spin-down rate and the inferred surface magnetic field, once a full timing solution is available. These measurements will shed light on the question of whether the fastest-spinning radio MSPs also typically have the lowest magnetic fields. The current limit of $B\lt 1.3\times {10}^{8}$ G already qualifies PSR J0952−0607 as one of the most weakly magnetized pulsars known.

Including PSR J0952−0607, there are 14 Galactic radio MSPs with $\nu \gt 500\,\mathrm{Hz}$, of which 5 are in black widow systems, 4 are isolated pulsars, 3 are in redback systems, and 2 have white dwarf companions. The abundance of black widow and redback systems among the fastest-spinning MSPs may hint at an evolutionary origin, possibly related to the accretion process and the amount of accreted matter (Hessels 2008; Papitto et al. 2014).²⁰ We note that radial velocity measurements and light curve modeling of the black widow companions to PSR B1957+20 ($\nu =622$ Hz) and PSR J1301+0833 ($\nu =543$ Hz) indicate that the MSPs are heavy, with masses of 2.40 ± 0.12 M${}_{\odot }$ (van Kerkwijk et al. 2011) and ${1.74}_{-0.17}^{+0.20}$ ${M}_{\odot }$ (Romani et al. 2016), respectively. Future photometric and spectroscopic observations of the companion of PSR J0952−0607 could test this hypothesis.

Besides the short spin period, PSR J0952−0607 is remarkable due to its steep radio spectrum. At $\alpha \sim -3$, its spectral index is among the steepest known compared to recent studies by Kuniyoshi et al. (2015), Kondratiev et al. (2016), and Frail et al. (2016). Our LOFAR and GBT observations should be robust against scintillation, as the bandwidths and integration times used substantially exceed the scintillation bandwidth and timescale toward PSR J0952−0607 (Cordes & Lazio 2002).

The steep spectrum of PSR J0952−0607 adds to the emergent picture where the fastest-spinning MSPs tend to have the steepest spectra, but also that the steepest spectra MSPs tend to be detected by Fermi in γ-rays (Kuniyoshi et al. 2015; Frail et al. 2016). Given that the fastest-spinning radio MSPs tend to have aligned radio and γ-ray profiles (Espinoza et al. 2013; Johnson et al. 2014), it is suggestive that these tendencies are pointing to a commonality in the radio and γ-ray emission mechanism, where the fast spin frequency leads to emission of γ-rays co-located with steep spectrum radio emission. Continued radio timing of PSR J0952−0607 will allow us to fold the Fermi γ-ray photons and test the alignment between the radio and γ-ray profiles.

The discovery of PSR J0952−0607 and the previous LOFAR discovery of PSR J1552+5437 (Bassa et al. 2017; Pleunis et al. 2017) demonstrate the potential of MSP searches at unconventionally low radio observing frequencies. Such searches are more sensitive to the population of ultra-steep-spectrum radio MSPs ($\alpha \lt -2.5$) and can explore whether even faster-spinning, potential submillisecond pulsars can be formed in nature but have previously been missed at higher observing frequencies. Unfortunately, low-frequency searches still suffer from limitations due to IISM scattering and eclipses. However, recent results have demonstrated the power of selecting such sources in low-frequency radio interferometric imaging surveys, where such effects do not affect the detectability of the source, and then following up with deep time-domain searches (Frail et al. 2016). Both imaging-led and direct time-domain search approaches should continue to be exploited.

While there is currently no clear theoretical expectation that the radio spectral index should depend on spin rate, we note that the shrinking size of the light cylinder ${r}_{c}=48(P/1\,\mathrm{ms})$ km could plausibly play a role in spectral steepening because the magnetospheric size becomes comparable to, or smaller, compared to the typical emission height. Furthermore, it is interesting to consider whether the radio emission from the fastest-spinning MSPs is dominated by giant pulse emission from a region co-located with the high-energy γ-ray emission.

We thank LOFAR Science Operations and Support for their help in scheduling and effectuating these observations. We also thank Caroline D'Angelo, Gemma Janssen, and Alessandro Patruno for useful discussions. C.G.B. and J.W.T.H. acknowledge support from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007–2013)/ERC grant agreement No. 337062 (DRAGNET; PI: Hessels). R.P.B. had received funding from the ERC under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 715051; Spiders). This Letter is based on data obtained with the International LOFAR Telescope (ILT) under projects LC7_002, DDT7_002 and LT5_003. LOFAR is the Low Frequency Array designed and constructed by ASTRON. It has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The Isaac Newton Telescope and its service programme are operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work was performed in part under DOE Contract DE-AC02-76SF00515.

Facilities: LOFAR - , ING:Newton - , Swift - , GBT - , Fermi. -

Software: cdmt, PRESTO, PSRCHIVE, Astropy, ESO-MIDAS, HEASOFT.