Detection and Timing of Gamma-Ray Pulsations from the 707 Hz Pulsar J0952−0607

The Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Atwood et al. 2009) has proven itself to be a powerful instrument in gamma-ray pulsar astronomy. Since its 2008 launch the LAT has been operating in an all-sky survey mode. LAT data are used to identify promising pulsar candidates for deep, targeted radio searches and find gamma-ray pulsations in blind or follow-up searches (for a review see, e.g., Caraveo 2014). The 10 year time span of the all-sky LAT data is also useful for establishing precise pulsar-timing ephemerides of new discoveries.

Radio pulsar searches targeting the sky positions of LAT sources have been very successful in finding isolated and binary millisecond pulsars (MSPs; e.g., Ray et al. 2012). The targeted sources are typically chosen to have three properties: (a) They are "unassociated," which means that the source has no plausible counterpart belonging to a known gamma-ray-emitting source class (e.g., Acero et al. 2015). (b) They have curved spectra. This is parametrized in the Fermi-LAT source catalogs by the curvature significance, determined by the difference in log-likelihood between spectral models with curved spectra (e.g., a log parabola or exponentially cutoff power law) versus power-law spectra (Nolan et al. 2012). For most gamma-ray pulsars, curved spectra are preferred with >95% confidence (e.g., Abdo et al. 2013). (c) They show only little variability in brightness over time, which is indicated in the Fermi LAT source catalogs by the variability index, the chi-squared of the monthly flux with respect to the average flux. In the Fermi LAT Third Source Catalog (3FGL; Acero et al. 2015), only 2 out of 136 pulsars had variability indices corresponding to significant variability above the 99% confidence level. Combined, the last two properties are good indicators for gamma-ray pulsars. However, we note that the transitional MSPs (for a review see, e.g., Jaodand et al. 2018) are an important exception, with significant changes in gamma-ray flux associated with transitions between accretion- and rotation-powered states (Stappers et al. 2014; Johnson et al. 2015).

Searches following this approach continue to find pulsars by using radio observing frequencies ν above 300 MHz. Pulsar surveys around 350 MHz are run by the Green Bank Telescope (GBT; Stovall et al. 2014) and the Arecibo telescope (Cromartie et al. 2016). The Giant Metrewave Radio Telescope searches around 607 MHz (Bhattacharyya et al. 2013). Another survey around 820 MHz is run by the GBT (Ransom et al. 2011). Finally Parkes (Camilo et al. 2015), Nançay (Cognard et al. 2011) and Effelsberg (Barr et al. 2013) search around 1.4 GHz. Radio observations at higher frequencies suffer less from dispersion (dispersion delay t_d ∝ ν⁻²) and scattering (scattering timescale τ_s ∝ ν^−4.4; Levin et al. 2016) but a pulsar's radio luminosity falls rapidly with observing frequency (radio flux density S_ν ∝ ν^α with spectral index −3.0 < α <−0.5 for most known pulsars; Frail et al. 2016a). At observing frequencies above 1.4 GHz scattering becomes negligible away from the Galactic Center and pulsars that are bright above this frequency can be useful for pulsar timing arrays (e.g., Verbiest et al. 2016; Tiburzi 2018).

However, there might be a population of steep-spectrum (α < −2.5) radio pulsars that are most easily detectable at frequencies below 300 MHz. Searches by Frail et al. (2018) for steep-spectrum sources within the localization regions of unidentified Fermi-LAT sources in continuum images from the Giant Metrewave Radio Telescope all-sky survey at 150 MHz led to the discovery of six new MSPs and one normal pulsar. These detections suggest that many steep-spectrum pulsars may have been missed by high-frequency radio surveys, which favor pulsars with flatter spectra (Bates et al. 2013). Additionally, some emission models suggest that pulsars' radio beams are wider at low frequencies (e.g., Story et al. 2007), making pulsars whose radio beams miss our line of sight at GHz frequencies potentially detectable at lower frequencies. Low-frequency radio observations of gamma-ray pulsars can therefore provide an additional test of the viewing-angle explanation for the large number of radio-quiet pulsars discovered by the LAT (e.g., Abdo et al. 2009; Wu et al. 2018). Indeed, one emission model for the recently discovered radio-quiet MSP PSR J1744−7619 (Clark et al. 2018) suggests that radio pulsations may only be detectable at low radio frequencies.

Pleunis et al. (2017) performed very-low-frequency pulsar searches at 115–155 MHz with the Low-Frequency Array (LOFAR; Stappers et al. 2011; van Haarlem et al. 2013). This was possible due to new semi-coherent de-dispersion techniques that mitigate the smearing due to dispersion (Bassa et al. 2017a). The searches targeted unassociated sources from the 3FGL catalog (Acero et al. 2015). An isolated MSP, PSR J1552+5437, was detected first in radio and subsequently in gamma-rays (Pleunis et al. 2017).

Bassa et al. (2017b) conducted another LOFAR survey using the same observing configuration. The 23 targets were unassociated gamma-ray sources selected from a Fermi-LAT source list constructed from 7 years of "Pass 8" LAT data (see Atwood et al. 2013).

In this survey they discovered PSR J0952−0607, a binary radio MSP with a spin frequency of 707 Hz (Bassa et al. 2017b). It is in a binary system with a very-low-mass companion star (M_c ∼ 0.02 M_⊙) with an orbital period of 6.42 hr. PSR J0952−0607 is the fastest-spinning known neutron star outside of a globular cluster: The only pulsar spinning faster (716 Hz) is PSR J1748−2446ad, which is located in the globular cluster Terzan 5 (Hessels et al. 2006). In contrast to pulsars in globular clusters, which experience significant but unknown acceleration due to the gravitational potential within the cluster (Prager et al. 2017), the intrinsic spin-down rate of PSR J0952−0607 can be measured directly. From this, pulsar properties like the dipole surface magnetic field strength and spin-down power can be inferred. These factors are thought to govern the poorly understood accretion and ablation processes through which binary systems containing a pulsar evolve (Chen et al. 2013). Measurements of the magnetic fields of rapidly spinning pulsars are important because the origin of the low magnetic field strength of MSPs is currently unexplained, with one popular theory being that the accreted matter buries the surface magnetic field. On the other hand recent work questions if this mechanism is effective enough (Mukherjee 2017).

To determine the pulsar properties requires precise timing solutions from frequent observations of a pulsar over several years. For some pulsar parameters (e.g., the spin frequency and spin-frequency derivative) the measurement uncertainty is directly related to the total span of observations. Furthermore, time spans shorter than 1 year cover less than a full cycle of the annual Roemer delay, introducing degeneracies between the spin frequency, spin-frequency derivative, and sky position. The radio-timing solution of PSR J0952−0607 reported by Bassa et al. (2017b) is based on observations spanning approximately 100 days, and thus suffers from these issues.

Radio searches targeting unassociated Fermi-LAT sources have been particularly successful at discovering "spider pulsars," a class of extreme binary pulsars with semi-degenerate companion stars (i.e., not neutron stars or white dwarfs). These systems are categorized as "black widows" if the companion star has extremely low mass (M_c ≪ 0.1 M_☉, as is the case for PSR J0952−0607) and as "redbacks" if the companion star is heavier (M_c ∼ 0.15–0.4 M_☉) (Roberts 2013). Optical light curves of these systems reveal that the pulsar emission heats the nearly Roche-lobe filling companion (Breton et al. 2013). Observations of orbitally modulated X-ray emission shows that interactions between the pulsar and companion star winds produce intra-binary shocks (e.g., Roberts et al. 2014).

For many spider pulsars the radio pulsations are completely absorbed by intra-binary material during parts of their orbit (e.g., Fruchter et al. 1988), indicating that the companion stars are also ablated by the pulsar. At low radio frequencies these eclipses can cover a large fraction of the orbit (e.g., Stappers et al. 1996; Archibald et al. 2009; Polzin et al. 2018), complicating radio-timing campaigns. In contrast, gamma-ray pulsations are essentially unaffected by eclipses.

A unique value of the LAT data is that a pulsar's discovery in gamma-rays often enables the immediate measurement of the pulsar parameters over the 10 year span in which the LAT has been operating. LAT data have been used to find precise timing solutions for many pulsars including radio-quiet and radio-faint pulsars (Ray et al. 2011; Kerr et al. 2015; Clark et al. 2017). In the case of PSR J2339−0533, a strongly eclipsing redback pulsar, gamma-ray timing was essential for building a coherent timing solution, and enabled the discovery of large variations of the orbital period (Pletsch & Clark 2015).

In this work we present the discovery and analysis of pulsed gamma-ray emission from PSR J0952−0607. The pulsar itself is very faint in gamma-rays, and required novel search and timing methods with greater sensitivity. The resulting timing ephemeris extends the rotational and orbital history of PSR J0952−0607 back 7 years to 2011. This allows us to determine the pulsar's spin-down power and surface magnetic field strength, making it the fastest known pulsar for which these measurements can be made.

The paper is organized as follows. In Section 2 we describe the pulsation search and detection within LAT data. The timing analysis and resulting timing solution for PSR J0952−0607 are presented in Section 3. New radio and optical observations as well as a search for continuous gravitational waves are discussed in Section 4. Finally, in Section 5 we discuss the implications of the results presented and we conclude in Section 6.

2.1. Data Preparation

2.2. Search

For many pulsars, LAT data covering several years of observation time are needed for significant pulsation detection (e.g., Hou et al. 2014). Searching for pulsations requires assigning every gamma-ray photon with the pulsar's rotational phase Φ (defined in rotations throughout the paper) at the time of emission. To do this a phase model Φ(t, λ) is used that depends on time t and (for circular-binary pulsars) on a set of at least seven parameters ${\boldsymbol{\lambda }}=(f,\dot{f},\alpha ,\delta ,{P}_{\mathrm{orb}},x,{t}_{\mathrm{asc}})$. These parameters are needed to (1) correct the photon arrival times for the LAT's movement with respect to the solar system barycenter (sky position α and δ), (2) in the case of a circular binary, account for the pulsar's movement around the center of mass (orbital period ${P}_{\mathrm{orb}}$, projected semimajor axis x, and epoch of ascending node ${t}_{\mathrm{asc}}$), and (3) describe the pulsar's rotation over time (spin frequency f and spin-frequency derivative $\dot{f}$).

The ephemeris obtained by timing a radio pulsar over a short interval ${T}_{\mathrm{obs}}$ often does not determine the parameters precisely enough to coherently fold the multiple years of LAT data. For ${T}_{\mathrm{obs}}\lt 1\,\mathrm{yr}$ the spin and position parameters of the pulsar are strongly correlated (i.e., degenerate). Over longer ${T}_{\mathrm{obs}}$ the uncertainties in the spin parameters scale with negative powers of ${T}_{\mathrm{obs}}$. The uncertainty in the orbital-period scales with ${T}_{\mathrm{obs}}^{-1}$ if ${T}_{\mathrm{obs}}\gg {P}_{\mathrm{orb}}$.

Searches for binary gamma-ray pulsars are therefore computationally expensive, as a multidimensional parameter space must be searched with a dense grid (Pletsch et al. 2012). The radio detection and timing are crucial to constrain the relevant parameter space that has to be searched to find the gamma-ray pulsations.

Using the radio data Bassa et al. (2017b) found that PSR J0952−0607 is in a circular-binary orbit. Furthermore, they measured α and δ by identifying the companion star using optical data taken with the Wide Field Camera (WFC) on the 2.5 m Isaac Newton Telescope on La Palma. Barycentering the radio data according to α and δ obtained from the optical data resulted in an upper limit on $\dot{f}$ and determined f more accurately. Furthermore the radio timing constrained the orbital parameters ${P}_{\mathrm{orb}}$, x, and ${t}_{\mathrm{asc}}$.

The gamma-ray pulsation search exploited preliminary constraints from radio timing of the pulsar combined with the optical position.

In the gamma-ray pulsation search we used the H statistic (de Jager et al. 1989). It combines the Fourier power from several harmonics incoherently by maximizing over the first M harmonics via

$$\begin{eqnarray}&&H=\mathop{\max }\limits_{1\leqslant M\leqslant {M}_{\max }}\left(4-4M+\displaystyle \sum _{n=1}^{M}{{ \mathcal P }}_{n}\right),\end{eqnarray} \tag{ 1 }$$

with M_max = 20 as suggested by de Jager et al. (1989). The Fourier power in the nth harmonic is given by

$$\begin{eqnarray}&&{{ \mathcal P }}_{n}=\displaystyle \frac{1}{{\kappa }^{2}}{\left|\displaystyle \sum _{j=1}^{N}{w}_{j}{e}^{-2\pi {in}{\rm{\Phi }}({t}_{j})}\right|}^{2},\end{eqnarray} \tag{ 2 }$$

with the normalization constant

$$\begin{eqnarray}&&{\kappa }^{2}=\displaystyle \frac{1}{2}\displaystyle \sum _{j=1}^{N}{w}_{j}^{2}.\end{eqnarray} \tag{ 3 }$$

The construction of a grid for this search was done using a distance "metric" on the parameter space (Balasubramanian et al. 1996; Owen 1996). This is a second-order Taylor approximation of the fractional loss in squared S/N due to an offset from the parameters of a given signal. The metric allows one to compute analytically the density of an optimally spaced grid. This method was successfully used in the blind search (i.e., a search for a previously undetected pulsar) for the black widow PSR J1311−3430 (Pletsch et al. 2012).

The metric components for the parameters of an isolated pulsar are given in Pletsch & Clark (2014), and the additional components required to search for a binary pulsar will be described in an upcoming paper (L. Nieder et al. 2019, in preparation). The grid point density computed with the metric varies throughout the parameter space. The grid density in α and δ increases as f increases. This is also the case for the orbital parameters. In addition, for ${P}_{\mathrm{orb}}$ and ${t}_{\mathrm{asc}}$ the grid point density increases with the projected semimajor axis, x. The small x typical for black-widow pulsars with their low-mass companions therefore greatly reduces the required density.

In addition, when performing a harmonic-summing search, any parameter offset results in a phase offset at the nth harmonic that is a factor of n larger than at the fundamental. To avoid this, the search grid density must be increased by a factor of M_max in each parameter. Fortunately, known gamma-ray pulsars have the most power in the first few harmonics (Pletsch & Clark 2014). We therefore designed the search grid to lose at most 1% of the Fourier power in the fifth harmonic in each dimension. The harmonic summing was also truncated at M_max = 5 to reduce computing cost. The required number of points in the search grid was reduced this way by a factor of 4⁵ (≈1000) compared to a grid built for M_max = 20. This search grid was designed to be very dense since the pulsar signal was expected to be weak due to the small number of photons.

Based on the distance metric we built a hypercubic grid covering the relevant parameter space in f, $\dot{f}$, α, δ, and ${P}_{\mathrm{orb}}$. This means that the parameter space is broken down into smaller cells. The edges of these cells are parallel to the parameter axes and of equal length in each dimension as computed by the distance metric. We note that a simple hypercubic grid is sufficient because the metric is nearly diagonal (off-diagonal terms are small; Nieder et al. 2019, in preparation), and the dimensionality is low. For higher dimensional parameter spaces hypercubic grids become extremely wasteful. The projected semimajor axis and the epoch of the ascending node were known precisely enough from the radio ephemeris that no search over these parameters was necessary. In summary, we performed a grid-based search over five parameters (f, $\dot{f}$, α, δ, and ${P}_{\mathrm{orb}}$), while keeping two parameters (x and ${t}_{\mathrm{asc}}$) fixed to the values from the radio-timing solution.

The search used 2 × 10⁵ CPU-core hours, meaning that the search would have taken 24 years to compute on a single core. Therefore, we distributed the work in chunks over  8000 CPU cores of the ATLAS computing cluster (Aulbert & Fehrmann2009), and the search took only 2 days.

2.3. Detection

3.1. Methods

We performed a timing analysis to measure precisely the parameters describing the pulsar's evolution over the observation time. We also allowed additional parameters to vary to test for measurable orbital eccentricity and proper motion of the binary. Instead of using a fixed search grid we use a Monte Carlo sampling algorithm to explore the parameter space around the signal parameters detected in the search. The general timing methods are also described by Clark et al. (2015, 2017), extending the methods developed by Ray et al. (2011) and Kerr et al. (2015). We enhanced these methods with the option to marginalize over the parameters of the template pulse profile as described in detail later in this section.

The starting point for the timing procedure is the construction of a template pulse profile, $\hat{g}({\rm{\Phi }})$, for which we used a combination of N_p symmetrical Gaussian peaks (Abdo et al. 2013)

$$\begin{eqnarray}&&\hat{g}({\rm{\Phi }})=\left(1-\displaystyle \sum _{i=1}^{{N}_{p}}{a}_{i}\right)+\displaystyle \sum _{i=1}^{{N}_{p}}{a}_{i}\,g({\rm{\Phi }},{\mu }_{i},{\sigma }_{i}).\end{eqnarray} \tag{ 4 }$$

The term $a\,g({\rm{\Phi }},\mu ,\sigma )$ denotes a wrapped Gaussian peak with amplitude a, peaked at phase μ with width σ:

$$\begin{eqnarray}&&g({\rm{\Phi }},\mu ,\sigma )=\displaystyle \frac{1}{\sigma \sqrt{2\pi }}\displaystyle \sum _{k=-\infty }^{\infty }\exp \left(-\displaystyle \frac{{\left({\rm{\Phi }}+k-\mu \right)}^{2}}{2{\sigma }^{2}}\right).\end{eqnarray} \tag{ 5 }$$

The phase at the first peak μ₁ is chosen to be the reference phase for the template. Phases of any other peak i are measured relative to the first peak as phase offset μ_i−μ₁ to avoid correlation with the overall phase. The template is fit to the weighted pulse profile obtained from the phase-folded data by maximizing over the likelihood

$$\begin{eqnarray}&&{ \mathcal L }(\hat{g},{\boldsymbol{\lambda }})=\displaystyle \prod _{j=1}^{N}[{w}_{j}\hat{g}({\rm{\Phi }}({t}_{j},{\boldsymbol{\lambda }}))+(1-{w}_{j})].\end{eqnarray} \tag{ 6 }$$

The Bayesian information criterion (BIC; Schwarz 1978) is used to choose the number of peaks by minimizing

$$\begin{eqnarray}&&\mathrm{BIC}=-2\mathrm{log}({ \mathcal L }(\hat{g},{\boldsymbol{\lambda }}))+k\mathrm{log}\left(\displaystyle \sum _{j=1}^{N}{w}_{j}\right),\end{eqnarray} \tag{ 7 }$$

where the number of free parameters in the model is denoted by k. Thus, adding a new parameter is penalized by $\mathrm{log}({\sum }_{j=1}^{N}{w}_{j})$ to avoid overfitting. The penalty factor for adding more Gaussian peaks to the template pulse profile scales with $k=3\times {N}_{{\rm{p}}}$ as each peak is described by three parameters.

As described by Clark et al. (2017), this template pulse profile is used to explore the multidimensional likelihood surface by varying the pulsar parameters with the goal to find the parameter combination that gives the maximum likelihood. We use our own implementation of the affine-invariant Monte Carlo method described by Goodman & Weare (2010) to run many Monte Carlo chains in parallel for the exploration and the efficient parallelization scheme described by Foreman-Mackey et al. (2013). The computations are distributed over several CPU cores.

This is repeated iteratively. Whenever a new best combination of parameters is found the template is updated using the new timing solution's phase-folded data. Usually this converges after a few iterations. Additional parameters (e.g., eccentricity) are added one after the other and the described timing procedure is restarted each time. Here again the BIC is used to decide whether the addition of a new parameter significantly improves the pulsar ephemeris. For the timing of bright pulsars (e.g., Clark et al. 2017) this iterative approach is sufficient.

For faint pulsars like PSR J0952−0607, the uncertainty in the gamma-ray pulse profile is not negligible. Using a fixed pulse profile template for weak pulsars could lead to systematic biases and underestimated uncertainties in the timing parameters. We therefore treated the template parameters in the same way as the pulsar parameters and let them vary jointly (as also done in An et al. 2017).

Joint variation of pulsar and template parameters results in larger but more realistic uncertainties on the pulsar parameters but should be used with a caveat. Varying pulsar parameters will always line up photons as close as possible to the same rotational phases to maximize the log-likelihood. The Monte Carlo algorithm finds combinations of parameters that lead to some photons being closer to the maximum of a peak and thus to a higher and narrower peak. But if these parameters do not describe the actual pulsar well, other photons will be shifted to phases outside the range of the peak, leading to a penalty preventing the acceptance of these parameter combinations. The joint variation of pulsar and template parameters however raises the chances of combinations that do not describe the actual pulsar well, as the peak position shifts to the phase where a combination of pulsar parameters leads to a narrow peak. This is a problem for a faint pulsar like PSR J0952−0607 as the penalty factor is weaker due to the smaller amount of photons. Furthermore for a pulsar like PSR J0952−0607 with two close peaks the penalty factor can be reduced by having one broader peak and one very narrow peak.

To address this problem we adjusted our priors on the template parameters. As for the pulsar parameters we used uniform priors for most template parameters. For the width parameters we used log-uniform priors and constrained them to peaks broader than 5% of a rotation, to disfavor extremely narrow peaks which only cover few photons, and narrower than half a rotation (full-width at half maximum ${\mathrm{FWHM}}_{i}\,=2\sqrt{2\mathrm{log}(2)}\,{\sigma }_{i}$ in the range 0.05 < FWHM_i < 0.5). This led to a steadier rise in H statistic over time and a pulse profile similar to what we get when folding the gamma-ray data with the updated radio-timing solution (see Section 4.1) reported in Table 2. In Figure 1 we show 100 pulse profile templates randomly picked from the resulting template parameter distribution.

Zoom In Zoom Out Reset image size

Table 2.  Properties of PSR J0952−0607 from Gamma-Ray and Radio Timing

Parameter	Gamma-Ray	Radio
Span of timing data (MJD)	55750a–58289	57759–58555
Reference epoch (MJD)	57980	57980
Timing Parameters
R.A., α (J2000.0)	09h52m08322(2)	09h52m0832141(5)
Decl., δ (J2000.0)	−06°07'2351(4)	−06°07'23490(2)
Spin frequency, f (Hz)	707.3144458307(7)	707.31444583103(6)
Spin-frequency derivative, ${\dot{f}}_{\mathrm{obs}}$ (Hz s−1)	−2.382(8) × 10−15	−2.388(4) × 10−15
Dispersion measure, DM (pc cm−3)		$22.411533(11)$
Orbital period, ${P}_{\mathrm{orb}}$ (day)	0.267461034(7)	0.2674610347(5)
Projected semimajor axis, x (lt-s)	0.0626670b	0.0626670(9)
Epoch of ascending node, ${t}_{\mathrm{asc}}$ (MJD)	57980.4479516b	57980.4479516(5)
Template Pulse Profile Parameters
Amplitude of first peak, α1	0.65(18)
Phase of first peak, μ1	0.431(39)
Width of first peak, σ1	0.064(23)
Amplitude of second peak, α2	0.35(24)
Phase offset of second to first peak, μ2−μ1	0.198(27)
Width of second peak, σ2	0.040(52)
Derived Properties (combined results)
Spin period, Pobs (ms)	1.414
Spin-period derivative,c ${\dot{P}}_{\mathrm{int}}$ (s s−1)	4.6 × 10−21
Characteristic age,d τc (Gyr)	4.9
Spin-down power,d $\dot{E}$ (erg s−1)	6.4 × 1034
Surface B-field,d ${B}_{\mathrm{surf}}$ (G)	8.2 × 107
Light-cylinder B-field,d BLC (G)	2.7 × 105
Galactic longitude, l (°)	243.65
Galactic latitude, b (°)	+35.38
NE2001 distance, (kpc)	${0.97}_{-0.53}^{+1.16}$
YMW16 distance, (kpc)	${1.74}_{-0.82}^{+1.57}$
Optical distance, (kpc)	${5.64}_{-0.91}^{+0.98}$
Gamma-ray luminosity,e Lγ (erg s−1)	3.1 × 1032 × (d/1 kpc)2

Notes. Numbers in parentheses are statistical 1σ uncertainties. The JPL DE405 solar system ephemeris has been used and times refer to TDB. Phase 0 is defined for a photon emitted at the pulsar system barycenter and arriving at the solar system barycenter at the reference epoch MJD 57,980.

^aValidity range of timing solution when the data starts at MJD 54,682. ^bFixed to values from radio-timing solution. ^cAssuming no proper motion, see Section 5. ^dProperties are derived as described in Abdo et al. (2013) on the basis of the estimated intrinsic spin-period derivative ${\dot{P}}_{\mathrm{int}}$. ^eAssuming no beaming and distance d = 1 kpc.

Download table as:  ASCIITypeset image

3.2. Solution

4.1. Updated Radio Timing

4.2. Optical Photometry

Bassa et al. (2017b) presented an r'-band light curve of the optical companion to PSR J0952−0607 taken by the WFC on the 2.5 m Isaac Newton Telescope on La Palma. The orbital light curve features a single maximum peaking at r' ≈ 22 at the pulsar's inferior conjunction, interpreted as being due to the pulsar heating the inside face (the "dayside") of a tidally locked companion. Bassa et al. (2017b) modeled this light curve with the Icarus package (Breton et al. 2012), finding that PSR J0952−0607 is likely to have an inclination angle i ∼ 40°, but the lack of color information precluded a robust estimate of other system parameters (e.g., companion temperature, heating, companion radius).

To more fully investigate the optical counterpart to PSR J0952−0607, we obtained multicolor photometry using ULTRACAM (Dhillon et al. 2007) on the 3.58 m New Technology Telescope (NTT) at ESO La Silla, and HiPERCAM (Dhillon et al. 2016, 2018) on the 10.4 m Gran Telescopio Canarias (GTC) on La Palma. The observation specifics are given in Table 3.

Table 3.  New Optical Photometry of the Companion of PSR J0952−0607

Night Beginning	Instrument+Telescope	Filters	${\phi }_{\mathrm{orb}}$	Airmass	Seeing	Photometric
2018 Jun 3	ULTRACAM+NTT	us, gs, is	0.64–1.09	1.1–2.1	10–20	yes
2018 Jun 4	ULTRACAM+NTT	us, gs, is	0.37–0.71	1.1–1.6	10–30	no
2019 Jan 12	HiPERCAM+GTC	us, gs, rs, is, zs	0.77–0.92	1.25–2.0	<15	yes
2019 Jan 13a	HiPERCAM+GTC	us, gs, rs, is, zs	0.37–0.72	1.25–2.0	15–20	no
2019 Mar 2b	ULTRACAM+NTT	us, gs, is	0.91–1.29	1.1–1.6	08–12	no
2019 Mar 3	ULTRACAM+NTT	us, gs, is	0.72–0.88	1.2–1.4	12–24	no
			1.16–1.72	1.1–1.9

Notes. Orbital phases are in fractions of an orbit, with ${\phi }_{\mathrm{orb}}=0$ corresponding to the pulsar's ascending node. The ULTRACAM data from 2018 were taken as a series of 20 s exposures in g_s and i_s, and 60 s in u_s. The 2019 ULTRACAM observations were taken with 10 s exposures in g_s and i_s, and 30 s in u_s. The HiPERCAM data cover u_s, g_s, r_s, i_s, and z_s simultaneously with exposure times of 60 s in u_s, g_s, r_s, and 30 s in i_s and z_s.

^aDuring an episode around ${\phi }_{\mathrm{orb}}=0.6$ seeing reached over 23 and 20 exposures had to be removed. ^bWe removed several frames due to intermittent clouds during the observations when the transmission dropped to nearly zero.

Download table as:  ASCIITypeset image

These data were calibrated and reduced using the ULTRACAM²⁰ and HiPERCAM²¹ software pipelines. Standard CCD calibration procedures were applied using bias and flat-field frames taken during each run.

We extracted instrumental magnitudes using aperture photometry, and performed "ensemble photometry" (Honeycutt1992) to correct for airmass effects and varying transparency. Magnitudes in ${g}_{s},{r}_{s},{i}_{s}$, and z_s²² were calibrated using comparison stars chosen from the Pan-STARRS1 (Chambers et al. 2016) catalog, after fitting for a color term accounting for differences between our filter sets and the Pan-STARRS1 filters. The HiPERCAM u_s observations were flux calibrated using zero-points derived from observations of two Sloan Digital Sky Survey (SDSS) standard stars (Smith et al. 2002) taken on 2019 January 11. The resulting HiPERCAM magnitudes for three nearby stars to PSR J0952−0607 were used to flux calibrate the ULTRACAM u_s data. Finally, the airmass- and ensemble-corrected count rates (C) were converted to AB flux densities according to our measured zero-point counts in each frame (C₀) by S_AB = 3631(C/C₀) Jy.

4.3. Optical Light-curve Modeling

As in Bassa et al. (2017b), the Icarus software was used to estimate parameters of the binary system. To do this, we performed a Bayesian parameter estimation using the nested sampling algorithm MultiNest (Feroz et al. 2013) via the Python package PyMultiNest (Buchner et al. 2014). Icarus produces model light curves by computing a grid of surface elements covering the companion star, and calculating and summing the projected line-of-sight flux from each element. Here the flux from each surface element was computed by integrating spectra from the Göttingen Spectral Library models of Husser et al. (2013).

In these fits we assumed that the companion star is tidally locked to the pulsar, and varied the following parameters: the companion star's "nightside" temperature (T_n); the "irradiating temperature" (T_irr defined such that the dayside temperature ${T}_{{\rm{d}}}^{4}={T}_{\mathrm{irr}}^{4}+{T}_{{\rm{n}}}^{4}$, under the assumption that the pulsar's irradiating flux is immediately thermalized and re-radiated, and therefore simply adds to the companion star's intrinsic flux at each point on the surface, as in Breton et al. 2013); the binary inclination angle (i); the Roche-lobe filling factor (f_RL, defined as the ratio between the companion's radius toward the pulsar and the inner Lagrange point (L1) radius); the distance modulus $(\mu =5{\mathrm{log}}_{10}(d)-5)$, with distance d in pc; and the mass of the pulsar (M_psr). At each point, the companion mass (M_c) and mass ratio (q = M_psr/M_c) were derived from the binary mass function according to the timing measurements of ${P}_{\mathrm{orb}}$ and x presented in Table 2. We also marginalize over interstellar extinction and reddening, parameterized by the E(B − V) of Green et al. (2018), scaled using the coefficients given therein for Pan-STARRS1 filter bands. We adopted a Gaussian prior for E(B − V) (truncated at zero), using the value from Green et al. (2018) for d > 1 kpc in the direction of PSR J0952−0607, E(B − V) = 0.065 ± 0.02, found by fitting the line-of-sight dust distribution using the apparent magnitudes of nearby main-sequence stars in the Pan-STARRS1 catalog. We adopted uniform priors on the remaining parameters (and uniform in $\cos i$), with M_psr and f_RL limited to lie within 1.2 < M_psr <2.5 M_⊙, and 0.1 < f_RL < 1. Temperatures T_n and T_d were constrained to lie within the range covered by the atmosphere models, 2300 < T < 12000 K.

At each point in the sampling, Icarus computed model light curves in each band. To account for remaining systematic uncertainties in the flux calibration, extinction, and atmosphere models, the model light curve in each band was re-scaled at each parameter location to maximize the penalized chi-squared log-likelihood. Overall calibration offsets were allowed for each band, and penalized by a zero-mean Gaussian prior on the scaling factor in each band with a width of 0.1 mag (a conservative estimate based on our calibration to the Pan-STARRS1 magnitudes). We also allowed small offsets between the calibrations for each ULTRACAM run and the HiPERCAM observations, which we penalized with an additional Gaussian prior with width 0.05 mag (also a conservative estimate from the differences in magnitudes of comparison stars in the field of view on each night). In initial fits, our best-fitting model resulted in a reduced chi-squared greater than unity. We therefore also re-scaled the uncertainties in each band to maximize the (re-normalized) log-likelihood at each point in the sampling. We also found that the fit improved substantially when we fit for a small orbital phase offset. Such orbital phase offsets are often seen in the optical light curves of black-widow pulsars and have been interpreted as being due to asymmetric heating from the pulsar, which could be caused by reprocessing of the pulsar wind by an intra-binary shock (e.g., Sanchez & Romani 2017).

The best-fitting light-curve model is shown in Figure 2, with posterior distributions for the fit parameters shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.4. Search for Continuous Gravitational Waves

The pulsar's spin period is defined as P = 1/f and the spin-period derivative is $\dot{P}=-\dot{f}/{f}^{2}$. The observed spin period for PSR J0952−0607 from gamma-ray and radio timing is P_obs = 1.414 ms and the observed spin-period derivative is ${\dot{P}}_{\mathrm{obs}}=4.76\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$.

The intrinsic spin-period derivative ${\dot{P}}_{\mathrm{int}}$ can be estimated from the observed value ${\dot{P}}_{\mathrm{obs}}={\dot{P}}_{\mathrm{int}}+{\dot{P}}_{\mathrm{Gal}}+{\dot{P}}_{\mathrm{Shk}}$. ${\dot{P}}_{\mathrm{Gal}}$ represents the part of the spin-period derivative caused by the relative Galactic acceleration (differential Galactic rotation and acceleration due to the Galactic gravitational potential; e.g., Damour & Taylor 1991; Nice & Taylor 1995), while ${\dot{P}}_{\mathrm{Shk}}$ accounts for the Shklovskii effect due to nonzero proper motion (Shklovskii 1970). Both contributions, ${\dot{P}}_{\mathrm{Gal}}$ and ${\dot{P}}_{\mathrm{Shk}}$, depend on the distance d to the pulsar.

The distance to PSR J0952−0607 is uncertain. The measured DM can be used to estimate the distance using Galactic electron-density models. The NE2001 model predicts ${0.97}_{-0.53}^{+1.16}\,\mathrm{kpc}$, while the YMW16 model predicts ${1.74}_{-0.82}^{+1.57}\,\mathrm{kpc}$. The uncertainties represent the 95% confidence regions (Yao et al. 2017). The model predictions of the DM as a function of d in the direction of the pulsar's sky position are shown in Figure 4. The models saturate at DM values that differ by ∼30% indicating the challenge and difficulty modeling the Galactic electron density. Still the distance predictions are consistent within the large uncertainty. On the other hand, the distance derived from optical modeling is ${5.64}_{-0.91}^{+0.98}\,\mathrm{kpc}$. This disagrees strongly with both DM distances and suggests that both DM models are overestimating the electron density in the direction of PSR J0952−0607. The distance discrepancy is discussed in more detail below.

Zoom In Zoom Out Reset image size

The estimated Galactic contribution is ${\dot{P}}_{\mathrm{Gal}}=(1.7,2.2,3.6)\,\times {10}^{-22}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ for the distance estimates d = (0.97, 1.74, 5.64) kpc. For the Shklovskii effect we then find the 95% confidence region to ${\dot{P}}_{\mathrm{Shk}}\in ([0,2.1],[0,3.8])\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ from the proper motion 95% confidence region (see Section 3.2) and for the (NE2001, YMW16) distances. The resulting 95% confidence region on ${\dot{P}}_{\mathrm{Shk}}$ for the optical distance exceeds past ${\dot{P}}_{\mathrm{obs}}$. Thus we only constrain the intrinsic spin-frequency derivative (at 95% confidence) to ${\dot{P}}_{\mathrm{int}}\in [2.44,4.59]\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ for the NE2001 model and ${\dot{P}}_{\mathrm{int}}\in [0.69,4.54]\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ for the YMW16 model. In the following, we conservatively assume zero proper motion (i.e., ${\dot{P}}_{\mathrm{Shk}}=0$) and used the fastest possible spin-down rate, ${\dot{P}}_{\mathrm{int}}=4.6\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$.

In Figure 5, PSR J0952−0607 is shown in a P–$\dot{P}$ diagram with the known pulsar population outside of globular clusters. The spin parameters of the more than 2000 radio pulsars are taken from the ATNF Pulsar Catalogue (see footnote 19) (Manchester et al. 2005).

Zoom In Zoom Out Reset image size

Furthermore we estimated the characteristic age τ_c, the spin-down power $\dot{E}$, the surface magnetic field strength ${B}_{\mathrm{surf}}$ and the magnetic field strength at the light cylinder B_LC (see Table 2). To calculate these values we assumed the pulsar to be a magnetic dipole with a canonical radius r_psr = 10 km and moment of inertia I_psr = 10⁴⁵ g cm² (e.g., Abdo et al. 2013). The same assumptions were used to plot the contour lines in Figure 5.

Despite spinning so rapidly, the gamma-ray energy flux of PSR J0952−0607 is on the fainter end of the gamma-ray MSP population. There are several reasons why gamma-ray pulsars might appear faint, including large distance, high background, or low luminosity (Hou et al. 2014). PSR J0952−0607 is not in a high-background region. The large distance derived from the optical modeling could be a possible explanation but disagrees with the distance estimates derived from the dispersion measure, d = (0.97, 1.74) kpc (NE2001, YMW16). The inferred gamma-ray luminosity is L_γ = 4πd²F_γ f_Ω ≈ 3.1 ×10³² × (d/1 kpc)² erg s⁻¹. The measured LAT energy flux F_γ is given in Table 1 and we assumed no beaming (i.e., f_Ω = 1). The gamma-ray efficiency is ${\eta }_{\gamma }={L}_{\gamma }/\dot{E}\approx 0.5 \% \times {(d/1\mathrm{kpc})}^{2}$. At the optical distance, η_γ ≈ 16% is typical of gamma-ray MSPs (Abdo et al. 2013), while at the DM-derived distance, η_γ ∼ 1% would be unusually low.

Due to the non-detection of PSR J0952−0607 in X-rays (F_X < 1.1 × 10⁻¹³ erg s⁻¹ cm⁻², Bassa et al. 2017b) we can only give a lower limit for the gamma-ray-to-X-ray flux ratio F_γ/F_X > 20. This limit is at the lower end of the observed distribution but still consistent with the literature (Marelli et al. 2011, 2015; Abdo et al. 2013; Salvetti et al. 2017).

The peak of the observed optical light curve is fairly broad in orbital phase. This requires either low inclination such that part of the heated face of the companion is visible over a large range of orbital phases, or for the companion to be close to filling its Roche lobe, such that the tidal deformation results in an "ellipsoidal" component peaking at ${\phi }_{\mathrm{orb}}=0.5$ and ${\phi }_{\mathrm{orb}}=1.0$ (with ${\phi }_{\mathrm{orb}}=0$ corresponding to the pulsar's ascending node) where the visible surface area of the companion is largest. Our best-fitting Icarus model favors the latter explanation, with f_RL ≈ 88% and i ≈ 61°. However, high filling factors imply a larger and hence more luminous companion, and therefore require greater distance, with our model having d ∼ 4.7–6.6 kpc.

We tried to refit the optical light curve with the distance fixed at the YMW16 distance of d = 1.74 kpc, but the resulting model has a significantly worse fit, and the low filling factor required results in an extremely high volume-averaged density for the companion (ρ) in excess of 100 g cm⁻³. For comparison, the densest known black-widow companions have densities of around 50 g cm⁻³ (e.g., PSR J0636+5128 Kaplan et al. 2018), with the record being that of the black-widow candidate 3FGL J1653.6−0158 in a 75 minute orbit (Romani et al. 2014) where ρ ≳ 70 g cm⁻³. These objects have been proposed to be the descendants of ultra-compact X-ray binaries, but this origin is unlikely for PSR J0952−0607 given its much longer orbital period (van Haaften et al. 2012). If the DM distances are assumed, the required density suggests that the companion star consists mostly of degenerate matter. A low filling factor may also explain the absence of radio eclipses seen from PSR J0952−0607. Alternatively, the low-density, large-distance solution has ρ ∼ 2.75 g cm⁻³, close to the density of brown dwarfs of similar mass and temperature given by the model considered in Kaplan et al. (2018).

We note that similar discrepancies in model distances were seen by Sanchez & Romani (2017) when using a direct-heating model. Romani & Sanchez (2016) and Sanchez & Romani (2017) considered models that additionally include a contribution from reprocessing of the pulsar wind by an intra-binary shock, which can wrap around the companion star. This can produce broader light curves for lower filling factors as some heating flux is redirected further around the sides of the companion star, and can also explain the small phase offset required for our direct-heating model by asymmetry in the shock front. Such a model may improve the fit for lower distances and filling factors, although an extremely high companion density would still be required to match the YMW16 distance. A likely explanation therefore could be that some heating flux is reprocessed by a shock, and the system has a moderate distance and filling factor, somewhat larger than required by the YMW16 value, but below those predicted by our direct-heating model. While more complex irradiation models (e.g., Romani & Sanchez 2016) may be required to address this issue, a full investigation of alternative models is beyond the scope of this study.

In both the small and large distance cases, we find that the nightside temperature of the companion is T_n ≈ 3000 ± 250 K at 95% confidence. We also find a well-constrained irradiating temperature of T_irr = 6100 ± 350 K, higher than that found from the single-band fit performed in Bassa et al. (2017b). This heating parameter can be compared to the total energy budget of the pulsar by calculating the "efficiency," ε, of conversion between spin-down power ($\dot{E}$) and heating flux (Breton et al. 2013)

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{4\pi {A}^{2}\sigma {T}_{\mathrm{irr}}^{4}}{\dot{E}},\end{eqnarray} \tag{ 8 }$$

with ε ∼ 20% being typical for black-widow systems. The efficiency is also shown in Figure 3, calculated from T_irr and from the orbital separation (A = x(1 + q)/sin i) at each point. We find that heating represents a larger fraction of the pulsar's total energy budget (ε ∼ 22% to 48% with 95% confidence) than the observed gamma-ray emission η_γ ≈ 0.5% × (d/1 kpc)². This estimate assumes that the pulsar's heating flux is emitted isotropically. As pointed out by Draghis & Romani (2018), some models of pulsar gamma-ray emission predict stronger beaming toward the pulsar's rotational equator, and an MSP's rotation should be aligned with the orbital plane as a result of the spin-up process. The actual gamma-ray luminosity directed toward the companion may therefore be higher than we observe. Our optical fits suggest a relatively face-on inclination (further evidenced by the lack of eclipses observed in radio observations, which often occur far outside the companion's Roche lobe), and so the comparative faintness of the pulsar's observed gamma-ray emission could be explained by the large viewing angle, and the fact that flux is preferentially emitted in the equatorial plane. A full modeling of the pulsar's phase-aligned radio and gamma-ray pulse profiles would provide an additional test of this scenario by estimating the viewing and magnetic inclination angles, and the relative beaming factors along our line of sight and in the equatorial plane. So far this is inhibited by the low significance of the gamma-ray light curve but with the continuing LAT mission this might be possible with more gamma-ray data in the future.

Alternatively, the difference between the heating flux and gamma-ray emission may suggest that another mechanism, e.g., the pulsar wind or intra-binary shock heating (Romani & Sanchez 2016; Wadiasingh et al. 2017), is responsible for heating the companion. Indeed, there is evidence for this being the case for the transitional PSR J1023−0038 where the optical heating is apparently unchanged between the MSP and low-mass X-ray binary (LMXB) states (Kennedy et al. 2018) despite a $5\times$ increase in the gamma-ray flux (Stappers et al. 2014).

As the optical counterpart to PSR J0952−0607 is faint (peaking at r' ≈ 22), it will be difficult to improve upon this picture of the system. While it may be possible to improve upon the dayside temperature measurement with optical spectroscopy in the future, the companion is effectively undetectable at minimum (r' > 25.0), precluding optical spectroscopic measurements of the companion's nightside temperature. We are also unable to constrain the mass of PSR J0952−0607 using the optical data. Constraining the pulsar mass would require a precise measurement of the binary mass ratio, which can be obtained for black-widow systems by comparing the radial velocities of the pulsar and companion. Unfortunately, the optical counterpart of PSR J0952−0607 is too faint (r' ∼ 23 at quadrature when the radial velocity is highest) for spectroscopic radial velocity measurements to be feasible even with 10 m class telescopes.

The gamma-ray source shows no significant variability as all flux measurements are consistent with the mean flux level. The calculated variability index also indicates a non-varying source. Here it is important to note that due to the low flux of the source the time bins had to be 750 days long to keep statistical precision. Therefore the variability index was calculated from only five independent time bins. Variations on shorter timescales can also not be found this way.

The gamma-ray pulse profile of PSR J0952−0607 shows two peaks that are separated by μ₂−μ₁ ≈ 0.2 rotations. This is typical for gamma-ray MSPs. More than half of them are double peaked with a peak separation of 0.2–0.5 rotations (Abdo et al. 2013). The radio pulse profile also shows two peaks with similar separation, with the radio pulse slightly leading the gamma-ray pulse (see Figure 6). The phase lag between the gamma-ray and radio pulse profile seems to be ∼0.15 (the majority of two-peaked MSPs show phase lags of 0.1–0.3; Abdo et al. 2013). Due to a covariance between f and dispersion measure (see Section 4.1) we were not able to measure significant variations in the dispersion measure. A change in dispersion measure of 10⁻³ pc cm⁻³ over the course of the Fermi mission would lead to an error in the phase offset of 13%.

Zoom In Zoom Out Reset image size

Gamma-ray pulsars are a good way to identify the maximum spin frequency of neutron stars. Among the 10 fastest Galactic field pulsars only 1 pulsar has not been detected in gamma-rays. Until the discovery of the 707 Hz pulsar PSR J0952−0607, the first MSP, PSR B1937+21, and the first black-widow pulsar, PSR B1957+20, were the fastest-spinning gamma-ray pulsars known (Guillemot et al. 2012). Still, the mass-shedding spin limit for neutron stars is typically placed much higher at around 1200 Hz (Cook et al. 1994; Lattimer & Prakash 2004). One mechanism that could prevent neutron stars from spinning up to higher frequencies is the emission of gravitational waves (for a recent work on this subject see, e.g., Gittins & Andersson 2019). Another option could be that the spin-up torque might be smaller for faster pulsars with lower magnetic field strengths (Patruno et al. 2012; Bonanno & Urpin 2015).

The estimated intrinsic spin-period derivative implies a very low surface magnetic field of 8.2 × 10⁷ G for PSR J0952−0607. Assuming nonzero proper motion would result in an even lower surface magnetic field estimate. Just nine pulsars, including the gamma-ray pulsar with the lowest surface magnetic field in the ATNF Pulsar Catalogue (see footnote 19) (Manchester et al. 2005), PSR J1544+4937 (Bhattacharyya et al. 2013), show lower inferred surface magnetic fields (Figure 7). The surface B-field of the other recent LOFAR-detected pulsar, PSR J1552+5437, is only slightly bigger (Pleunis et al. 2017). This might be a hint that pulsars with low B-fields also have steeper radio spectra.

Zoom In Zoom Out Reset image size

The pulsar distribution in Figure 7 indicates a lower limit on the magnetic field strength independent of the spin frequency. The equilibrium spin period as predicted by Alpar et al. (1982) is ${P}_{\mathrm{eq}}\propto {B}_{\mathrm{surf}}^{6/7}\,{R}_{\mathrm{psr}}^{18/7}\,{M}_{\mathrm{psr}}^{-5/7}\,{\dot{M}}_{\mathrm{accr}}^{-3/7}$ with pulsar radius R_psr, mass M_psr, and accretion rate ${\dot{M}}_{\mathrm{accr}}$, which indicates that the lowest spin periods can be reached for low magnetic field strengths and high accretion rates. Nevertheless high accretion rates lead to a rapid decrease of the magnetic field strength and for low magnetic field strengths the angular momentum transfer is slower (Bonanno & Urpin 2015). In order to spin up to millisecond periods a limiting magnetic field strength and accretion rate can be set as a result of the amount of time a neutron star can spend accreting matter being limited by the age of the universe (Pan et al. 2018). For a neutron star with a mass of 1.4 M_⊙, a radius of 10 km and a minimum accretion rate of 7.26 × 10⁻¹¹ M_⊙ yr⁻¹ we get a minimum magnetic field strength of B_surf ≳ 3.3 × 10⁷ G, which is consistent with the observed pulsar population.

No continuous gravitational waves are detected from PSR J0952−0607, which is to date the fastest-spinning pulsar targeted for gravitational wave emission. The 95% upper limit on the intrinsic gravitational wave amplitude is set to ${h}_{0}^{95 \% }=6.6\times {10}^{-26}$. The corresponding upper limit on the ellipticity is ${\epsilon }^{95 \% }=3.1\,\times {10}^{-8}\times (d/1\,\mathrm{kpc})\times ({10}^{45}\,\,{\rm{g}}\,{\mathrm{cm}}^{2}/I)$, where I is the principal moment of inertia of the pulsar. The intrinsic gravitational wave amplitude at the detector needed to account for all of the spin-down energy lost due to gravitational wave emission is ${h}_{0}^{\mathrm{sd}}\,=1.5\times {10}^{-27}\times {(1\mathrm{kpc}/d)\times (I/{10}^{45}{\rm{g}}{\mathrm{cm}}^{2})}^{1/2}$, corresponding to an ellipticity of ${\epsilon }^{\mathrm{sd}}=7.0\times {10}^{-10}\times (1\,\mathrm{kpc}/d)$.

As for many other high-frequency pulsars, the indirect spin-down upper limit on h₀ is smaller and more constraining than our measured gravitational wave upper limit, in this case by a factor of ≈45 at 1 kpc. For a more likely larger distance the factor would be even greater, so it is not surprising that a signal was not detected (Abbott et al. 2019). The quoted spin-down upper limit could be inaccurate if the measured spin down were affected by radial motions, if the distance were smaller than estimated or if the moment of inertia of the pulsar were different than the fiducial value of 10⁴⁵ g cm². In the case of PSR J0952−0607 it is unlikely that all these effects could bridge a gap of nearly two orders of magnitude, but in line with the "eyes-wide-open" spirit of previous searches for gravitational waves from known pulsars (see Abbott et al. 2019, 2017; Aasi et al. 2014 and references therein) we all the same perform the search.

Using a sensitive, fully coherent pulsation search technique, we detected gamma-ray pulsations from the radio pulsar PSR J0952−0607 in a search around the parameters reported by Bassa et al. (2017b). New timing methods were developed to cope with the low signal strength, allowing us to measure the spin rate, sky position, and orbital period with high precision, and in agreement with the updated radio-timing ephemeris. Furthermore thanks to the longer gamma-ray time span we reliably constrained the intrinsic spin-period derivative ${\dot{P}}_{\mathrm{int}}\lesssim 4.6\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$. This measurement provides estimates of physical parameters such as the spin-down luminosity ($\dot{E}\lesssim 6.4\times {10}^{34}$ erg s⁻¹), and a surface magnetic field (${B}_{\mathrm{surf}}\lesssim 8.2\times {10}^{7}$ G) among the lowest of any detected gamma-ray pulsar. Although the resulting timing solution spans 7 years to the present data, we were unable to extend this to cover data earlier than MJD 55,750. We investigated several possible reasons. Flux variations could lead to the loss of pulsations. A time-varying orbital period as seen in several spider pulsars would cause a loss of phase coherence. With our current data we are not able to ascertain the true reason. In the absence of orbital-period variations or state changes, improved timing precision from additional data should help determine the cause.

We also obtained new multiband photometry of the pulsar's optical counterpart, and modeled the resulting light curve. To explain the observed optical flux, our models require either a much larger distance (∼5 kpc) than the DM-distance estimates of 0.97 kpc (NE2001) to 1.74 kpc (YMW16), or a small and extremely dense companion ρ ≫ 100 g cm⁻³. More complex optical models including intra-binary shocks might help to solve this discrepancy, but a full investigation of other models is beyond the scope of this work. We found that the pulsar flux heating the companion star accounts for a much larger fraction of the pulsar's spin-down power (∼50%) than is converted to observed gamma-ray emission (0.5% at 1 kpc), although this difference is reduced if our larger distance estimate is adopted.

Despite the extensive analysis of PSR J0952−0607 and its companion, the study of this pulsar has not ended as some questions remain unanswered. The LAT and LOFAR continue to take gamma-ray and radio data on this source, and we plan to obtain more optical data.

LAT gamma-ray data has helped to find many new MSPs by providing promising candidates (Ray et al. 2012). Sophisticated methods to identify more pulsar candidates within LAT sources have been developed (e.g., Lee et al. 2012; Saz Parkinson et al. 2016). For instance, Frail et al. (2016b) identified 11 promising MSP candidates by checking for steep-spectrum radio sources coincident with LAT sources. With the approach successfully used in this paper, new binary MSP candidates can be searched for pulsations and upon detection the pulsar can be precisely timed within months after its discovery. Identifying more of the rapidly rotating spider pulsars will be helpful to study further the observed neutron star parameter limits like the maximum spin frequency and the minimum surface magnetic field strength.

We thank the referee for pointing out the uncertainties in the DM/distance models, and suggesting the arguments given around Figure 4. This work was supported by the Max-Planck-Gesellschaft (MPG) and the ATLAS cluster computing team at AEI Hannover. C.J.C., R.P.B., and D.M.-S. acknowledge support from the ERC under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 715051; Spiders). C.G.B. and J.W.T.H. acknowledge support from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 337062 (DRAGNET; PI: Hessels). This work was supported by an STSM Grant from COST Action CA16214. M.R.K. is funded through a Newton International Fellowship provided by the Royal Society. Work at NRL is supported by NASA.

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. This work performed in part under DOE Contract DE-AC02-76SF00515.

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.

Part of this work is based on data obtained with the international LOFAR Telescope (ILT) under project codes LC7_018, DDT7_002, LT5_003, LC9_041, and LT10_004. LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage 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 ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d'Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK.

HiPERCAM and V.S.D. are funded by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) under ERC-2013-ADG grant agreement number 340040 (HiPERCAM). ULTRACAM and V.S.D. are funded by the UK Science and Technology Facilities Council. This work is based on observations made with the GTC, installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, on the island of La Palma, and on observations made with ESO Telescopes at the La Silla Paranal Observatory.

This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.

Software: Fermi Science Tools, dspsr (van Straten & Bailes 2011), psrchive (Hotan et al. 2004), tempo2 (Hobbs et al. 2006; Edwards et al. 2006), NE2001 (Cordes & Lazio 2002), YMW16 (Yao et al. 2017), MultiNest (Feroz et al. 2013), PyMultiNest (Buchner et al. 2014), ULTRACAM/HiPERCAM software pipelines, Icarus (Breton et al. 2012), psrqpy (Manchester et al. 2005; Pitkin 2018), Astropy (Astropy Collaboration et al. 2013, 2018), matplotlib (Hunter 2007), NumPy (Oliphant 2006; van der Walt et al. 2011).

It is important to estimate the false-alarm probability P_FA to know if the gamma-ray detection is real. As described in Section 2.3, there is no known analytical expression for the false-alarm probability of the maximum value from an H statistic search over a dense, multidimensional parameter grid. Deriving the probability distribution for the maximum value of a multidimensional "random field" is difficult and approximate solutions are only known for simple cases such as Gaussian or chi-squared random fields (Adler & Taylor 2007). While the power in a single harmonic does follow a chi-squared random field in the presence of random noise, the known solutions cannot be applied in this case due to the maximization over summed harmonics and penalty factors defining the H statistic, and the fact that the metric density varies between different summed harmonics. Even for chi-squared random fields, there is no simple "trials factor" that can be applied to the single-trial false-alarm probability (which for the H statistic was derived by Kerr 2011): the false-alarm probability depends on the volume, shape, and dimensionality of the search space (Adler & Taylor 2007). A full discussion of this is beyond the scope of this work. Below, we show empirically that a simple trials factor approach overestimates the detection significance, and describe the "bootstrapping" method that we used to overcome this.

The false-alarm probability for a single H statistic trial is

$$\begin{eqnarray}&&{P}_{\mathrm{FA}}({H}_{{\rm{m}}}\,| \,a)={e}^{-a{H}_{{\rm{m}}}},\end{eqnarray} \tag{ 9 }$$

with scaling factor a = 0.3984 (de Jager & Büsching 2010; Kerr 2011). This formula can be used to estimate the significance of the maximum H statistic value after n independent trials

$$\begin{eqnarray}&&{P}_{\mathrm{FA}}({H}_{{\rm{m}}}\,| \,a,n)=1-{\left[1-{e}^{-a{H}_{{\rm{m}}}}\right]}^{n}.\end{eqnarray} \tag{ 10 }$$

We assume at first that our search contained a number of "effective" independent trials (N_eff) that is some unknown fraction of the number of actual trials (i.e., the number of grid points at which we evaluated the H statistic). We then estimated N_eff from the results of our search as follows. We divided our parameter space into n_seg = 2 × 17 × 13 = 442 segments in f, $\dot{f}$, and ${P}_{\mathrm{orb}}$ respectively. The number of segments in f and $\dot{f}$ is determined by the parameter space volumes, which were searched in parallel, as only the highest H statistic values from each were stored. To ensure that all segments were independent from the pulsar signal, we removed all grid points within those segments which were close (according to the parameter space metric; see Section 2.2) to the pulsar parameters.

The highest H statistic of each of the segments is plotted in the normalized histogram in Figure 8. We fit for the effective number of trials (as done by, e.g., Kruger et al. 2002) by maximizing the likelihood,

$$\begin{eqnarray}&&L(n,a| {H}_{m,i})=\displaystyle \prod _{i}p({H}_{m,i}| a,n)\end{eqnarray} \tag{ 11 }$$

for our set of H statistic values, according to the probability density function for H_m after n trials (the derivative of Equation (10))

$$\begin{eqnarray}&&p({H}_{{\rm{m}}}| a,n)=a\,n{\left[1-{e}^{-a{H}_{{\rm{m}}}}\right]}^{n-1}\,\exp (-a\,{H}_{{\rm{m}}}).\end{eqnarray} \tag{ 12 }$$

However, as shown in Figure 8, the tail of the best-fitting distribution is significantly underestimated, leading to overestimated significances for large H statistic values. This demonstrates that there is no simple effective trials factor that can be applied to estimate the overall significance.

Zoom In Zoom Out Reset image size

To overcome this, we performed a second fit, maximizing over the likelihood for both n and a. The resulting best-fitting distribution is also shown in Figure 8. We found the best-fitting scaling factor to be $\hat{a}\approx 0.284$, meaning the probability density function is flatter and gives a more conservative estimate for the significance. We note that this should not apply in general, and will depend, among other factors, on the dimensionality of the search space and the number of harmonics summed.

Finally, we use $\hat{a}$ and multiply the best-fitting n (the best-fitting per-segment trials factor) by n_seg, and apply Equation (10) to obtain an approximation to the false-alarm probability for the maximum H statistic value. For the candidate pulsar signal, this was P_FA = 0.33%. For comparison the candidate with the largest H statistic from a segment of the search not affected by the pulsar signal had P_FA = 56%.