UvA-DARE (Digital Academic Repository) Detection and Timing of Gamma-Ray Pulsations from the 707 Hz Pulsar J0952-0607

The Low-Frequency Array radio telescope discovered the 707 Hz binary millisecond pulsar ( MSP ) J0952 − 0607 in a targeted radio pulsation search of an unidenti ﬁ ed Fermi gamma-ray source. This source shows a weak energy ﬂ ux of F γ  =  2.6  ×  10 − 12 erg cm − 2 s − 1 in the energy range between 100 MeV and 100 GeV. Here we report the detection of pulsed gamma-ray emission from PSR J0952 − 0607 in a very sensitive gamma-ray pulsation search. The pulsar ’ s rotational, binary, and astrometric properties are measured over 7 years of Fermi -Large Area Telescope data. For this we take into account the uncertainty on the shape of the gamma-ray pulse pro ﬁ le. We present an updated radio-timing solution now spanning more than 2 years and show results from optical modeling of the black-widow-type companion based on new multiband photometric data taken with HiPERCAM on the Gran Telescopio Canarias on La Palma and ULTRACAM on the New Technology Telescope at ESO La Silla ( based on observations collected at the European Southern Observatory, Chile; programme 0101.D-0925, PI: Clark, C. J. ) . PSR J0952 − 0607 is now the fastest-spinning pulsar for which the intrinsic spin-down rate has been reliably constrained ( ) . The inferred surface magnetic ﬁ eld strength of is among the 10 lowest of all known pulsars. This discovery is another example of an extremely fast spinning black-widow pulsar hiding within an unidenti ﬁ ed Fermi gamma-ray source. In the future such systems might help to pin down the maximum spin frequency and the minimum surface magnetic ﬁ eld strength of MSPs.


Introduction
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-rayemitting 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 chisquared 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 gammaray 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 , 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 ∝ν −2 ) 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 steepspectrum 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 viewingangle 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 radioquiet MSP PSR J1744−7619 ) 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 gammarays (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 e ) 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 semidegenerate 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 . 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 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 .
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 2 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.

Data Preparation
The gamma-ray source targeted by Bassa et al. (2017b) resulting in the detection of the radio pulsar PSR J0952−0607 and its optical counterpart (R.A. α J2000.0 = 09 h 52 m 08 319, decl. δ J2000.0 = −06°07′23 49) was discovered using 7 years of LAT data, but was too faint to be included in the 3FGL catalog (i.e., in 4 years of data; Acero et al. 2015). It is included in the successive 4FGL catalog based on 8 years of data as 4FGL J0952.1−0607 (The Fermi-LAT collaboration 2019).
To search for gamma-ray pulsations from PSR J0952−0607, we used "Pass 8" (Atwood et al. 2013) LAT data recorded between 2008 August 4 and 2017 January 19, consisting of SOURCE-class photons above 500 MeV instead of the standard 100 MeV. Since the LATʼs angular resolution for photons improves with energy (∼3.6 times higher angular resolution at 500 MeV compared to 100 MeV), we conservatively used 500 MeV to avoid potential contamination by other nearby sources not included in the 3FGL catalog. 14 The photons were selected using gtselect from the Fermi Science Tools 15 if they were within 10°of the celestial position of the optical counterpart to PSR J0952−0607, with a maximum zenith angle of 90°. Photons were only used if the LAT was in nominal science mode and if the rocking angle was below 52°. After these cuts 114706 LAT photons remained for further analysis. The analysis was performed using the P8R2_SOURCE_V6 instrument response functions (IRFs).
The sensitivity of a pulsation search can be greatly improved by weighting the contribution of each photon by its probability of having originated from the candidate pulsar (Bickel et al. 2008;Kerr 2011). The weights are computed based on the LAT response function and a spectral model of a point source. They are used in the search and the timing analysis for background suppression without the need for arbitrary position or stronger energy cuts.
To produce the necessary spectral model we performed a binned spectral analysis with gtlike. We added a putative pulsar source with an exponentially cutoff power law to represent its spectrum (Nolan et al. 2012) fixed to the position of the pulsar's optical counterpart reported by Bassa et al. (2017b). We used the templates gll_iem_v06.fits for the Galactic diffuse emission (Acero et al. 2016) and iso_P8R2_SOURCE_V6_v06.txt 16 for the isotropic diffuse background. The spectral analysis included all 3FGL sources within 15°of the pulsar position and the spectral parameters for point sources within 5°of the target were allowed to vary.
For each photon within 5°of the pulsar's optical position a probability weight w j was calculated with gtsrcprob. To reduce the computing cost of the search, we only included photons with w j >3.1%. This weight cutoff value was chosen such that only 1% of the expected pulsation signal-to-noise ratio (S/N) would be lost. After applying the cutoff N=1354 actual or å = w 193.7 j "effective" photons remain. Upon the detection of PSR J0952−0607, we performed a dedicated spectral analysis with an extended data set in order to enhance the pulsation significance and to model its spectral characteristics more precisely. We used the same event selection and IRFs (see above) but accepted photons without cuts on the rocking angle as this cut was found to be overly conservative. 17 We extended the data set to include photons between 2008 August 4 and 2018 June 21. We lowered the threshold of photon energies down to 100 MeV to further constrain the spectral characteristics. We used the Preliminary LAT 8-year Point Source List 18 (FL8Y) to construct our source model. The FL8Y source associated with the pulsar, FL8Y J0952.2−0608, was replaced by a point source fixed to the position of the detected gamma-ray pulsar. All FL8Y sources within 15°of the pulsar position were included and the spectral parameters for point sources within 5°of the pulsar were allowed to vary.
We computed the residual TS map to search for noncataloged weak gamma-ray sources in the vicinity of the pulsar. The test statistic = -  TS 2 log source log no source ( ( ) ( )) quantifies how significant a source emerges from the background, where the likelihood  of a model with and without a source is compared (Nolan et al. 2012;Acero et al. 2015). Six uncatalogued sources with TS>10 (∼3σ) within 5°of the pulsar position were found and added to the source model. Using this new source model we reran the analysis. The result of the spectral analysis for PSR J0952−0607 is shown in Table 1. Here, we also give TS cut which is computed like TS but comparing an exponentially cutoff power-law model and a power-law model without cutoff (Abdo et al. 2013).
In the timing analysis we used all photons with weights w j >1.5%, which is chosen as in the search such that 99% of Note.Gamma-ray spectrum based on LAT data between MJD 54,682-58,289 over the standard energy range from 100 MeV to 100 GeV. The first reported uncertainties are statistical, while the second uncertainties are systematic, determined by re-analyzing the data with bracketing IRFs and artificially changing the normalization of the Galactic diffuse model by ±6%, as described in Abdo et al. (2013).
the S/N remains. This leaves N=4642 actual or å = w 331.4 j effective photons.

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 l a d = f f P x t , , , , , , orb 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 orb , projected semimajor axis x, and epoch of ascending node t asc ), and (3) describe the pulsar's rotation over time (spin frequency f and spin-frequency derivative f˙).
The ephemeris obtained by timing a radio pulsar over a short interval T obs often does not determine the parameters precisely enough to coherently fold the multiple years of LAT data. For < T 1 yr obs the spin and position parameters of the pulsar are strongly correlated (i.e., degenerate). Over longer T obs the uncertainties in the spin parameters scale with negative powers of T obs . The uncertainty in the orbital-period scales with -T obs 1 if T P obs 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 f˙and determined f more accurately. Furthermore the radio timing constrained the orbital parameters P orb , x, and t 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 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 orb and t 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 5 (≈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, f˙, α, δ, and P 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, f˙, α, δ, and P orb ), while keeping two parameters (x and t asc ) fixed to the values from the radio-timing solution.
The search used 2×10 5 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 & Fehrmann 2009), and the search took only 2 days.

Detection
To ensure that the signal was inside the covered parameter space we searched over wide ranges in the highly correlated f (4σ), α, and δ (5σ each), where σ is the parameter uncertainty obtained from preliminary radio and optical observations. The chosen search range for P orb (3σ) was smaller because the radiotiming-derived P orb was not degenerate with the other parameters.
Surprisingly, the largest H statistic appeared close to the edge of our search range in f and with a significant offset in α and δ. The latter was determined to be due to an error in the initial astrometric calibration of the optical images of the optical counterpart. After the discovery of this error only the corrected α and δ values were published by Bassa et al. (2017b). The offset in f arose from the strong correlation with α and δ. Therefore we started another search with the same settings starting from the highest f covered in the first search. The largest H statistic was H m =86.7 (without refining the parameters any further) and lay well within the combined search parameter space.
While this H statistic was far larger than any other found in our search, it is not easy to estimate the statistical significance (or false-alarm probability) of the maximum value found in a dense, multidimensional H statistic search (see Appendix). We therefore applied a "bootstrapping" procedure (described in the Appendix) to estimate the detection significance from the search results themselves, finding a trials-corrected false-alarm probability of P FA ≈3.3×10 −3 . After extending our data set to cover the extra year of data as explained in Section 2.1, and without using a weight cut (which is only introduced for computational reasons), we found that the H statistic value increased to H=102.9 without further refinement (i.e., in a single trial). Since no additional trials have been performed in this step, we can multiply our false-alarm probability estimate by the known single-trial false-alarm probability (Kerr 2011) for , giving an overall false-alarm probability of P FA ≈5.3×10 −6 in the extended data set, confirming the detection.

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. (2015Clark et al. ( , 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, F ĝ ( ), for which we used a combination of N p symmetrical Gaussian peaks (Abdo et al. 2013) å å The term m s F a g , , ( )denotes a wrapped Gaussian peak with amplitude a, peaked at phase μ with width σ: The phase at the first peak μ 1 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 −μ 1 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 The Bayesian information criterion (BIC; Schwarz 1978) is used to choose the number of peaks by minimizing where the number of free parameters in the model is denoted by k. Thus, adding a new parameter is penalized by ) to avoid overfitting. The penalty factor for adding more Gaussian peaks to the template pulse profile scales with =ḱ N 3 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 = FWHM i s 2 2 log 2 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.

Solution
Our timing solution is shown in Table 2. We did not find clear pulsations in the beginning of the Fermi mission at MJD 54,682 and therefore our timing solution starts at MJD 55,750 (see Figure 1). We discuss the absence of pulsations prior to MJD 55,750 below.
The gamma-ray pulse profile is likely double peaked as the double-peaked template is favored by the BIC over the singlepeaked template. The template parameters leading to the highest likelihood are given in Table 2.
All of the measured parameters are consistent with the initial published radio solution. The published values and uncertainties on α and δ from the optical counterpart are consistent and comparable to the ones in the gamma-ray timing solution (Bassa et al. 2017b). As expected from the much longer timing baseline the uncertainties on f and P orb are much smaller than in the initial radio-timing solution. Furthermore, it is possible to measure the spin-frequency derivative, 10 Hz s 15 1 . A second spin-frequency derivative, f¨, is clearly disfavored by the BIC. The gamma-ray timing solution is consistent with an updated radio ephemeris based on radio data spanning 796 days, and the parameter uncertainties are comparable or smaller (see Section 4.1 and Table 2).
It is not possible for us to confidently determine the proper motion as we find hints for and against nonzero proper motion. Allowing proper motion to vary jointly with the template parameters results in a significantly improved H statistic, loglikelihood, and BIC. The timing analysis sets the 95% confidence region on proper motion to m d Î -- 2 is 14.8 mas yr −1 with a 95% upper limit of 25.3 mas yr −1 . Typically, however, it is assumed that the H statistic rises linearly with exposure time and nonzero proper motion resulting from this timing analysis leads to a bumpier rise in the H statistic over time. This indicates that the proper motion resulting from our analysis might not be correct. Keeping the template fixed to the template parameters found by folding the gamma-ray data with the radio ephemeris results in a 95% confidence region on proper motion consistent with zero. Zero proper motion is also favored by the BIC. The same is found when using a single-peaked profile in the timing analysis and varying the template parameters jointly. 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. a Validity range of timing solution when the data starts at MJD 54,682. b Fixed to values from radio-timing solution. c Assuming no proper motion, see Section 5. d Properties are derived as described in Abdo et al. (2013) on the basis of the estimated intrinsic spin-period derivative P inṫ . e Assuming no beaming and distance d=1 kpc.
The upper limit on proper motion corresponds to a transverse velocity of v t =μ t d=120 km s −1 ×(d/1 kpc). This results in high, but not unrealistic transverse velocities when using the distances inferred from the dispersion measure (d = 0.97 kpc (Cordes & Lazio 2002 or d=1.74 kpc (Yao et al. 2017, hereafter YMW16)). As 90% of the known MSPs in the Australia Telescope National Facility (ATNF) Pulsar Catalogue 19  show transverse velocities below 200 km s −1 the proper motion upper limit is unrealistic for the higher distances predicted by the optical observations (4.7-6.6 kpc; see Section 5).
Unsurprisingly, we were unable to detect a significant timing parallax. The maximum parallax time delay for the abovementioned distance estimates is D »p t d 500 lt s 2 ,max 2 ( -) ( ) m 1 s. In comparison the resolution with which we can measure the arrival time of the pulse is Δμ 1 /f≈61 μs.
A circular orbit is clearly favored over an eccentric orbit by the BIC. The 95% upper limit on eccentricity is set to e<0.004.
The missing pulsations before MJD 55,750 seem odd as the tracks are clearly visible later in the mission (Figure 1). As the pulsar is not very bright one explanation might be Poisson variations in the flux leading to the loss of pulsations for a few hundred days. Possible pulsations before this period might be too weak to be picked up again as the phase uncertainty grows quickly outside the timing span. At the start of the mission (MJD 54,682) the phase uncertainty is ∼0.6 rotations, which could be a plausible explanation for loss of coherence.
In order to understand the nature of the non-detection of gamma-ray pulsations before MJD 55,750, we measured the gamma-ray flux of PSR J0952−0607 over time by sliding a 750 day long window in steps of 50 days over the LAT data. In each of these steps we calculated the gamma-ray flux of PSR J0952−0607 over the 750 days width of the window, which allowed us to measure the spectral parameters with reasonable precision. We found that the flux of PSR J0952 −0607 is lower in the beginning of the Fermi mission but the lower fluxes agree with the flux uncertainties from the full time span. The TS values follow the same trend as the gamma-ray fluxes in the sliding windows.
The gamma-ray source is too faint to test it unambiguously for variability. The windows need to cover 750 days to keep statistical precision. But that leaves only five independent time bins to calculate the variability index with Equation(4) from Abdo et al. (2010). The variability index computed with these five bins is 7.18 with 4 degrees of freedom, which is below the 99% confidence level of 13.277.
We also checked if the smaller 35°rocking angle used during the first year of the Fermi mission decreases the pulsation significance. However, the small rocking angle is actually favorable as the exposure for PSR J0952−0607 is ∼20% higher in the beginning of the mission.
Variations of the orbital period might be another reasonable explanation for the loss of clear pulsations. Such orbital-period variations have been measured for several spider pulsars, e.g., for the original black-widow pulsar PSR B1957+20 (Arzoumanian et al. 1994). Nevertheless the penalty for adding orbital-period derivatives led to an increase in the BIC. Similarly, no significant semimajor-axis derivative was found.

Updated Radio Timing
Observations of PSR J0952−0607 with LOFAR have been ongoing using an identical observational setup as in Bassa et al.  The phase-connected timing solution from Bassa et al. (2017b) was improved by using all LOFAR HBA observations that used 78 MHz of bandwidth (hence excluding the discovery and initial follow-up observations which used half the bandwidth). Pulse time-of-arrival (TOA) measurements were obtained by referencing pulse profiles of eight frequency channels per observation to a single analytic pulse profile template. This procedure presumes that our data are not sensitive to pulse profile shape variations with frequency, which was double-checked through inspection of the difference profiles of the top and bottom parts of the bandpass: no significant structures were detected. The analytic pulse profile was created using the PSRCHIVE (van Straten et al. 2012) package PAAS and was constructed from five von Mises functions that were fitted to the integrated body of observations and fully modeled any detectable pulse shapes. The resulting timing solution extends the timing baseline to 2.2 years and breaks the degeneracy between the astrometric and rotational parameters (see Table 2). Upon inspection of the data, a new covariance was detected, namely, between a significant (>4σ) decrease in the dispersion measure of this pulsar (which was found to be decreasing by 5 × 10 −5 pc cm −3 yr −1 ) and the spin period. Notwithstanding the significance of this decrease, the strong anticorrelation of this parameter with the pulse period suggests an underestimate of its measurement significance, which is commonly found in pulsar-timing analyses (e.g., Coles et al. 2011), particularly in nonperiodic parameters such as linear gradients in dispersion measure. Consequently this decrease was not included in our present analysis, but future monitoring to allow more robust disentanglement of the spin period and the dispersion measure variability is warranted. We find no evidence for radio eclipses in the six LOFAR observations with orbital phases between 0.15<f orb <0.35. Using the TOAs from this orbital phase range we set a 3σ upper limit on time delays due to additional dispersion of Δt<2.3 μs, and hence ΔDM<1. 2×10 −5 pc cm −3 . 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 , 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).

Optical Photometry
To more fully investigate the optical counterpart to PSR J0952−0607, we obtained multicolor photometry using ULTRACAM  Table 3.
These data were calibrated and reduced using the ULTRA-CAM 20 and HiPERCAM 21 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" (Honeycutt 1992) to correct for airmass effects and varying transparency. Magnitudes in g r i , , s s s , and z s 22 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 zeropoint counts in each frame (C 0 ) by S AB =3631(C/C 0 ) Jy.

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 ). 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  , 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 m = d 5 log 5 10 ( ( ) ), 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 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-ofsight 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 i cos ), with M psr and f RL limited to lie within 1.2<M psr < 2.5 M e , and 0.1<f RL <1. Temperatures T n and T d were 20 http://deneb.astro.warwick.ac.uk/phsaap/software/ultracam/html/ 21 http://deneb.astro.warwick.ac.uk/phsaap/hipercam/docs/html/ 22 ULTRACAM and HiPERCAM use identical higher-throughput versions of the SDSS filter set, which we refer to as Super-SDSS filters: u s , g s , r s , i s , and z s (Dhillon et al. 2018). 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.

Search for Continuous Gravitational Waves
We carried out a search for near-monochromatic continuous gravitational waves phase locked at twice the pulsar rotation phase for the source PSR J0952−0607 using data from the first and second runs (O1 23 and O2 24 ) of the two Advanced LIGO detectors (Vallisneri et al. 2015). The observation period spans 707 days from 2015 September to 2017 August and comprises 183 days (169 days) of data from the Hanford (Livingston) detector.
We employ the coherent multi-detector detection statistic  2 (Jaranowski et al. 1998;Cutler & Schutz 2005) that we  Notes.Orbital phases are in fractions of an orbit, with f = 0 orb 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 . a During an episode around f = 0.6 orb seeing reached over 2 3 and 20 exposures had to be removed. b We removed several frames due to intermittent clouds during the observations when the transmission dropped to nearly zero.
implemented in the LIGO-LALSUITE library. 25  2 is the log-likelihood maximized over the amplitude parameters i y h , cos , 0 and Φ 0 for a near-monochromatic 26 gravitational wave signal with given frequency and frequency-derivative values, from a source in a binary at a given sky position and with given orbital parameters, in Gaussian noise. h 0 is the intrinsic gravitational wave amplitude at the detector, ι the angle between the total angular momentum of the pulsar and the line of sight to it from Earth, ψ is the gravitational wave polarization angle and Φ 0 the signal phase at a nominal reference time. In this search we assume the gravitational wave frequency and frequency derivatives equal to twice the values measured for the pulsar rotation frequency and its derivatives. In Gaussian noise the detection statistic  2 follows a χ 2 -distribution with 4 degrees of freedom and non-centrality parameter equal to 0: the expected value is μ=4.0, and the standard deviation is s = 2 2. If a signal is present, the noncentrality parameter is proportional to the square of the intrinsic gravitational wave amplitude at the detector, h 0 , and to the total observation time.
The search yields the value =  2 9.9, which is well within the bulk of the distribution consistent with a null result. Based on the measured value of the detection statistic, we set a frequentist 95% upper limit on the intrinsic gravitational wave amplitude, h 0 95% , following a now standard procedure first developed by some of us (Abbott et al. 2004). h 0 95% is the smallest intrinsic gravitational wave amplitude such that 95% of the population of signals that could be emitted by PSR J0952 −0607 27 would yield a detection statistic value greater than the measured one, =  2 9.9. We find =´h 6.6 10 0 95% 26 . The uncertainty on this upper limit is~14%, including instrument calibration errors (Cahillane et al. 2017).

Discussion
The pulsar's spin period is defined as P=1/f and the spinperiod derivative is = -P 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 =´--P 4.76 10 s s obs 21 1 . The intrinsic spin-period derivative P inṫ can be estimated from the observed value = + + P P P P obs int Gal Shk˙˙˙. P 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 P Shk accounts for the Shklovskii effect due to nonzero proper motion (Shklovskii 1970). Both contributions, P Gal and P 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 electrondensity models. The NE2001 model predicts . 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.
The estimated Galactic contribution is =Ṕ . In Figure 5, PSR J0952−0607 is shown in a 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) .
Furthermore we estimated the characteristic age τ c , the spindown power Ė, the surface magnetic field strength B 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 45 g cm 2 (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 2 F γ f Ω ≈3.1× 10 32 ×(d/1 kpc) 2 erg s −1 . 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 h = »ǵ g L E d 0.5% 1 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 −13 erg s −1 cm −2 , 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  (Yao et al. 2017). To illustrate the discrepancy with these distance predictions, the 95% confidence region from the optical modeling is shown. The vertical, dashed-dotted line indicates the distance favored by the optical modeling. 27 The possible signals span uniformly distributed values of i -  1 cos 1 and of 0ψ2π. 11 distribution but still consistent with the literature (Marelli et al. 2011(Marelli et al. , 2015Abdo 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 f = 0.5 orb and f = 1.0 orb (with f = 0 orb 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 −3 . For comparison, the densest known black-widow companions have densities of around 50 g cm −3 (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 −3 . 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 −3 , 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 (Ė) and heating flux 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) 2 . 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 Figure 5. Spin period P and spin-period derivative Ṗ of the known pulsar population outside of globular clusters. The inset shows a zoomed-in view of the known MSP population. Isolated radio pulsars (light-gray pluses), binary radio pulsars (dark-gray squares), isolated gamma-ray pulsars (light-green crosses) and binary gamma-ray pulsars (dark-green circles) are shown. The subject of this paper, the gamma-ray pulsar PSR J0952−0607, is marked by an orange star. The lines denote constant characteristic age τ c (dotted), spin-down power Ė (dashed) and surface magnetic field strength B surf (dashed-dotted).
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 lowmass X-ray binary (LMXB) states (Kennedy et al. 2018) despite a5 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 μ 2 −μ 1 ≈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 −3 pc cm −3 over the course of the Fermi mission would lead to an error in the phase offset of 13%.
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 . 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 7 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) , PSR J1544+4937 (Bhattacharyya et al. 2013), show lower inferred surface magnetic fields (Figure 7). The surface B-field of the other recent LOFARdetected 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.
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 µ with pulsar radius R psr , mass M psr , and accretion rate M accṙ , 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 Figure 6. Aligned integrated gamma-ray and radio pulse profiles of PSR J0952 −0607 over two identical rotations. The black curve shows the weighted LAT photon counts after MJD 55,750 in a histogram with 30 bins per rotation. The green error bars show the phase uncertainty of the gamma-ray pulse profile. The estimated background level is indicated by the dashed blue line. The radio profile as seen by the LOFAR telescope in a 78 MHz band centered at 149 MHz is drawn in red. The error bars drawn in dark red indicate the possible phase shift of the radio pulse profile due to a dispersion measure variation of 10 −3 pc cm −3 over the time span of the Fermi mission. of 1.4 M e , a radius of 10 km and a minimum accretion rate of 7.26×10 −11 M e yr −1 we get a minimum magnetic field strength of B surf 3.3×10 7 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 6.6 10 0 95% 26 ( ) . As for many other high-frequency pulsars, the indirect spindown upper limit on h 0 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 45 g cm 2 . 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-wideopen" spirit of previous searches for gravitational waves from known pulsars (see Abbott et al. 2019Abbott et al. , 2017Aasi et al. 2014 and references therein) we all the same perform the search.

Conclusions
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 derivativé . This measurement provides estimates of physical parameters such as the spin-down luminosity ( E 6.4 10 34 erg s −1 ), and a surface magnetic field ( B 8.2 10 surf 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 −3 . 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  European Union's Seventh Framework Programme (FP7/ 2007(FP7/ -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: 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 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, f˙, and P orb respectively. The number of segments in f and 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.
15 Figure 8. Normalized histogram showing the highest H statistics for 442 subsets of our search space after excluding results affected by the pulsar signal. The dotted green and dashed blue curves show normalized probability density functions for the maximum H statistic obtained after n effective trials. The curves gave maximum likelihood after varying over n with fixed single-trial scaling factor a=0.3984 (dotted green) and after varying a and n jointly (dashed blue). The maximum H statistic for the pulsar H m =86.7 is marked by the vertical orange line. The red line (dashed-dotted) shows the false-alarm probability depending on H m computed with Equation (10) with a and n from the joint variation.