Discovery, Timing, and Multiwavelength Observations of the Black Widow Millisecond Pulsar PSR J1555$-$2908

One of the great successes of Fermi LAT mission (Atwood et al., 2009) has been the revelation that a large number of point sources in the GeV gamma-ray sky are powered by pulsars, including large populations both of energetic young pulsars and older, faster spinning millisecond pulsars. In the early days of the Fermi mission, these discoveries came via two techniques: (1) folding the observed gamma-rays at the pulse period of known radio pulsars (e.g. Abdo et al., 2009a) and (2) direct searches for periodicities in the gamma-ray photon arrival times themselves (e.g. Abdo et al., 2009b). Each of these techniques had its own biases and selection effects that left many gamma-ray pulsars undiscovered. The next highly successful technique was to use the locations of LAT gamma-ray sources that had pulsar-like characteristics (e.g. curved spectra and low variability) as targets for directed radio pulsar searches, performed by the Fermi Pulsar Search Consortium (PSC) (Ray et al., 2012). This led to a bounty of millisecond pulsar discoveries that has continued as the LAT survey observations reveal ever-more gamma-ray sources (Ransom et al., 2011; Cognard et al., 2011; Kerr et al., 2012; Ray et al., 2013; Bhattacharyya et al., 2013; Camilo et al., 2015; Cromartie et al., 2016). Most recently, a new technique was employed to reveal pulsar candidates among the LAT sources: Frail et al. (2016) and Frail et al. (2018) exploited the nearly full-sky 150 MHz TIFR GMRT Sky Survey (TGSS) radio survey performed by the Giant Metrewave Radio Telescope (GMRT) to identify steep spectrum radio point sources within the error regions of LAT sources. These candidates were quickly followed up by the PSC with deep radio searches, leading to several new millisecond pulsar (MSP) discoveries (initially reported in Frail et al. 2018) and demonstrating the power of this technique. Here we provide details on the discovery of PSR J1555$-$2908 and follow up observations that determined its orbit, discovered gamma-ray pulsations, searched for X-ray pulsations, and made initial studies of its optical companion star.

PSR J1555$-$2908  was initially identified from an image-based search for potential pulsar candidates. The method is described in more detail in Frail et al. (2018) but here we give a short summary. We searched for compact, steep spectrum radio sources within the error ellipses of unassociated Fermi sources from a preliminary version of the LAT 8-year catalog (hereafter 4FGL; Abdollahi et al., 2020). We used the GMRT 150 MHz All-Sky Radio Survey (TGSS ADR1; Intema et al., 2017) and the 1.4 GHz NRAO VLA Sky Survey (NVSS; Condon et al., 1998) to calculate initial two-point spectral indices and source angular diameters for all radio sources within Fermi unassociated sources. Outside of the NVSS survey declination limit ($-45^{\circ}$), but north of the declination limit of TGSS ($-53^{\circ}$), we used the 843 MHz Sydney University Molonglo Sky Survey (SUMSS; Bock et al., 1999) for the higher frequency of the two-point spectral index. For those rare compact, steep spectrum radio sources identified via this method, follow-up interferometric observations at arcsecond resolution were carried out in order to eliminate false positives from the sample, such as high redshift radio galaxies.

This image-based method differs from past searches of Fermi unassociated sources (e.g. Kerr et al., 2012) in that pulsar candidates are selected without regard to the properties of the Fermi source, such as its spectral shape, temporal variability, the size of the error ellipse, or its sky distribution. A single unresolved radio source was identified within the error ellipse of P86Y3595 (now known as 4FGL J1555.7$-$2908) with a spectral index $\alpha=-2.5\pm$0.2 (where $\alpha$ is given as $S_{\nu}\propto\nu^{\alpha}$). The radio source position was measured to be R.A.=$15^{\rm h}55^{\rm m}40^{\rm s}\!.69$, decl.=$-29^{\circ}08^{\prime}29^{\prime\prime}\!.0$ (J2000) with an uncertainty of about 2″.

As part of a long term effort to discover radio pulsars associated with Fermi sources, organized by the Fermi PSC, we made targeted observations of the candidates identified in Frail et al. (2018). We observed P86Y3595 on 2017 February 1 (MJD 57785) for 30 minutes at 820 MHz with the Robert C. Byrd Green Bank Telescope (GBT). We used the GUPPI pulsar backend to record 200 MHz of bandwidth with 2048 channels at 61.44 $\mu$s resolution. We performed an acceleration search of these data over a range of trial dispersion measures using PRESTO (Ransom et al., 2002) and identified a strong candidate pulsation with 1.79 ms period (559.4 Hz) at DM 75.91 pc cm^-3. Following the discovery, a 5 minute GBT observation at S-band (on MJD 57806) confirmed the discovery and determined the pulse width at that frequency to be a very narrow 3%.

We also obtained a 30-minute observation on MJD 57788 using the Karl G. Jansky Very Large Array (VLA) in phased-array pulsar mode. The data were taken while the array was being reconfigured; 20 antennas were in their D-configuration positions (maximum baseline 1 km), and 7 still in A-configuration (maximum baseline 30 km). The observation was split into three 10-minute scans in order to allow for re-phasing. Interferometric phases were determined using the nearby calibrator J1554$-$2704, then individual antenna data streams were coherently summed to create a high time resolution single-pixel data stream at the previously-determined radio source position. The summed voltage data were processed in real time using dspsr (van Straten & Bailes, 2011) and recorded using 2048 0.5-MHz channels (1024 MHz total bandwidth, centered on 1500 MHz) at 32 $\mu$s time resolution. It is worth noting that this VLA observation is quite broadband and exhibits profile evolution across the band. This profile evolution likely biases the measured DM from that observation.

We also detected the pulsation in a 2-minute GBT observation at 350 MHz that had been taken on MJD 56907 as part of the GBNCC survey (Stovall et al., 2014). Lastly we obtained a 5 hour GBT observation on MJD 57833 revealing an eclipse with a sharp ingress and clear pulse delays for about 10 minutes after egress. The eclipse duration is approximately 10% of the duration of the orbit, and the likely dispersive pulse delays are typical for black-widow pulsar systems. This observation was taken in coherently-dedispersed search mode, providing a calibrated flux density and high time resolution. A fault in the GUPPI backend at the time of the observation prevented us from extracting reliable polarization measurements from this observation. The eclipse and high S/N total intensity pulse profile from that observation are shown in Figure 2.

We also obtained five observations of varying duration using the 76-m Lovell telescope between 2019 November 15 and 2020 September 21. We used a backend based on ROACH FPGA boards (Bassa et al., 2016a) to record coherently de-dispersed, phase-folded data covering a bandwidth of 384 MHz in 0.25 MHz channels centred at 1534 MHz, with 256 phase bins per pulse period. To obtain approximate flux-density calibrations for these observations, we estimated the system equivalent flux density from observations of the Crab nebula on nearby epochs, adjusting for ground spillover due to the low elevation. We estimate the fractional uncertainty on this calibration to be around 20%. The pulsed flux density was then estimated from each observation by fitting the high S/N 820MHz GBT pulse profile as a template, along with an arbitrary baseline, to the flux-calibrated profiles. The average flux density over the five observations was $S_{1534}=0.20\pm 0.05$ mJy. A log of all the radio observations is presented in Table 1. The pulse profile at a range of frequencies is shown in Figure 1.

Combining these radio data, we obtained an initial orbital solution with a period of 5.60 hours and a semimajor axis of 0.151 lt-s. Because we have only a small amount of radio data, most of the system parameters are better determined by the gamma-ray timing in Section 4.2. However, the precise radio TOAs over the 5 hour observation allow us to determine some of the orbital parameters more precisely than is possible from the gamma-ray data alone. In particular, we measure an eccentricity $e=4.5\pm 1.5\times 10^{-6}$, $x=a_{1}\sin i=0.1514468(1)$ lt-s, and an epoch of ascending node of $t_{\text{asc}}=57785.53936388(3)$ MJD(TDB).

The initial (not phase-connected) radio timing parameters provided the seed necessary to make a computationally-tractable pulsation search in the gamma-ray data, as we describe here.

Although the distance inferred from the dispersion measure is large (see Section 7), PSR J1555$-$2908 is one of the 10 fastest MSPs in the Galactic field (i.e., outside of a globular cluster), and the LAT timing indicates that the spin down power is very high, at $3.1\times 10^{35}$ erg s^-1. This makes it a good candidate to search for X-ray pulsations, as other high $\dot{E}$ MSPs (e.g. B1937+21, B1821$-$24, and J0218+4232) exhibit bright non-thermal pulsations.

Motivated by this, we made a 122 ks observation of this source with the NICER X-ray telescope (NICER proposal #2527), with data accumulated from 2019 May 26 (ObsID 2527010101) through 2019 September 9 (ObsID 2527010142). Ray et al. (2019) provide a description of the NICER X-ray Timing Instrument and details of the data reduction for millisecond pulsar pulsation searches. We follow a similar procedure here, but since these data were taken with a large angle between the target and the Sun, optical loading was not a significant problem, so we did not mask any of the 52 active detectors. We analyzed our data with HEASoft version 6.27.2 (NICER tools version 2020-04-23_V007a) and updated the gain calibration to version nixtiflightpi20170601v006. Our initial data extraction included energies 0.25–10.0 keV, and made standard data cuts to exclude the SAA and ensure NICER was tracking the source and the source was $>20^{\circ}$ above the Earth limb. We made no cuts on magnetic cutoff rigidity, or the total count rate, in our initial extraction. This initial extraction yielded 117.8 ks of good time. For each photon, we computed the pulse phase using photonphase from PINT (Luo et al., 2021) and the timing model provided in Table 2.

To exclude high background regions, instead of making an arbitrary cut on count rate, we developed a tool ni_Htest_sortgti.py that divides the data into segments of no more than 100 s (and no less than 10 s). These segments are sorted by mean count rate. Since the count rate from the pulsar is very low and presumably constant, this is equivalent to sorting by background rate. Then the $H$-test detection statistic is evaluated cumulatively, going from the lowest background to the highest. This algorithm thus finds the optimal background rate cut that maximizes the S/N of any detected signal. The script repeats this process over a large grid of $E_{\mathrm{min}}$ and $E_{\mathrm{max}}$ energy cuts to search for an optimal energy band as well. For our detection statistic, we chose to use the $H$-test (de Jager et al., 1989), which is preferred to the $Z^{2}$ test for unknown pulse profiles that may be very sharp, as is seen in the non-thermal X-ray MSPs.

For PSR J1555$-$2908, this procedure did not reveal any very strong candidate pulsation. The highest $H$-test found corresponded to a single-trial significance of 3.05 $\sigma$, when searching the range 0.26–2.74 keV, with the algorithm selecting the lowest background 106.7 ks out of the total. In those data, the mean count rate was 0.80 s^-1. Since there were a large number of trials over the energy and background cuts, this does not represent a significant detection.

In the 0.25–2.0 keV band, the median count rate (source+background) was 0.75 s^-1, with 90% of segments below 1.13 s^-1. In the 2–8 keV band, the median count rate (source+background) was 0.23 s^-1, with 90% of segments $<0.6$ s^-1.

We performed differential time-series photometry at the position of PSR J1555$-$2908 using the 2.4-m Hiltner telescope of the MDM Observatory on 2018 June 12 and 15. A back-illuminated SITe CCD with $1024\times 1024$ $24\,\mu$ pixels, each subtending $0.\!^{\prime\prime}275$, was exposed through an SDSS $r$ filter. Integrations of 5 minutes each were obtained for a total of 7.16 hours over the two nights. The differential photometry was calibrated using a comparison star from Pan-STARRS. A variable star with a maximum brightness of $r=20.4$ and displaying a typical heating-dominated light curve was detected at the image-based radio position (Section 2), also consistent with the Gaia counterpart position and magnitude (Section 4.2) and the gamma-ray timing position (Table 2). A finding chart is shown in Figure 4.

It was immediately apparent that the timescale and pattern of variability are consistent with orbital modulation of a low-mass companion star heated by the pulsar wind. The phased light curve using the orbital ephemeris of Table 2 is shown in Figure 5. As expected, the peak brightness occurs near phase 0.75, superior conjunction of the companion, where we are viewing the heated side of the star. The “night” side of the companion was not detected; only upper limits of $r>23$ were obtained at phases $\phi=0.25\pm 0.1$.

The radio discovery of PSR J1555$-$2908 continues the trend of successfully finding radio pulsations from MSPs by targeting LAT sources. It also demonstrates the power of using imaging radio surveys to discover steep-spectrum point sources associated with LAT sources. These are prime pulsar candidates.

The pulsar discovery provides some constraints on the distance to the source. The dispersion measure of 75.9 pc cm^-3 gives distance estimates of 2.65 kpc for the NE2001 electron density model (Cordes & Lazio, 2002) and 7.55 kpc for the YMW16 model (Yao et al., 2017). The difference of the two estimates can be traced to the properties of the thick disk component in the electron density models. The DM at high Galactic latitude is dominated by the thick disk, which in the YMW16 model has a smaller central density and a larger scale height than in the NE2001 model, and an integrated column density that is slightly more than half that of NE2001. High-latitude pulsars are thus assigned a larger DM distance in YMW16. PSR J1555$-$2908 is at Galactic coordinates $(\ell,b)=(344.\!^{\circ}5,+18.\!^{\circ}5)$, and has the largest DM of neighboring pulsars within a $5^{\circ}$ radius. This pushes it to a large, albeit uncertain distance in YMW16.

From our flux-calibrated observations with the Lovell telescope, we estimated an L-band pulsed flux density of $S_{1534}=0.20\pm 0.05$ mJy. Comparing this to the imaging flux density measured at 150 MHz by Frail et al. (2018) of $S_{150}=133.9\pm 14.4$ mJy, and assuming a simple power-law spectrum $S_{\nu}\propto\nu^{\alpha}$, yields a two-point spectral index estimate of $\alpha=-2.8\pm 0.1$, confirming the steep spectrum.

PSR J1555$-$2908 is a fast and energetic pulsar, currently the 13th fastest galactic MSP in the ATNF Pulsar Catalogue. The small upper limit on the measured proper motion of 6.4 mas yr^-1 sets an upper limit on the Shklovskii correction to the spindown rate of at most 3%. The larger of the two DM distance estimates gives an intrinsic spindown luminosity of $2.98\times 10^{35}$ erg s^-1. This can be compared with the measured gamma-ray luminosity. The associated LAT source in the 4FGL DR3 catalog is 4FGL J1555.7$-$2908 with an energy flux above 100 MeV (G100) of $(4.66\pm 0.61)\times 10^{-12}$ erg cm^-2 s^-1. This corresponds to an efficiency of converting spindown energy into $>100$ MeV gamma-rays of 1.3% for $d=2.65$ kpc and 11% for $d=7.55$ kpc. These are both plausible values for gamma-ray MSPs (see Abdo et al., 2013). Alternatively, Abdo et al. (2013) define a rough heuristic relationship $L_{\gamma}^{h}=\sqrt{10^{33}\dot{E}}$ erg s^-1, which gives an expected gamma-ray luminosity of $1.8\times 10^{34}$ erg s^-1. This implies a distance of 5.7 kpc, but with substantial uncertainty based on the observed scatter around that relationship.

The X-ray pulsation upper limits can be compared to the population of gamma-ray millisecond pulsars whose properties are compiled in Abdo et al. (2013). The ratio of $>100$ MeV gamma-ray flux ($G_{100}$) to non-thermal X-ray flux ($F_{X}^{\mathrm{nt}}$; unabsorbed 0.3–10 keV) ranges over at least 2 orders of magnitude, from $\sim 20$ to over several thousand. If PSR J1555$-$2908 is like the energetic MSP PSR B1937+21 (essentially the most optimistic case), which has a $G_{100}/F_{X}^{nt}=23$, then $F_{X}$ could be as high as $3\times 10^{-13}$ erg s^-1. Assuming $n_{\rm H}=2.2\times 10^{21}$ cm^-2 scaled from the DM per He et al. (2013) and a power-law index of 2.0, yields a NICER count rate of 0.044 s^-1. The detectability of pulsations depends on the unknown pulse shape (ranging from the worst case of a pure sinusoid to a very narrow, low duty cycle pulse), but we can estimate the detection significance using the method of Buccheri et al. (1987). For this flux, we should have detected these pulsations if the source pulsed fraction was larger than 75% (for a sinusoid) or 40% (for a narrow pulse). For less favorable gamma-ray to X-ray flux ratios, e.g. with predicted count rates $<0.018$ s^-1, our observation would not have detected pulsations even with very narrow pulse shapes. Consequently, the X-ray observation is not sensitive enough to place very stringent constraints on the X-ray luminosity, but could have detected pulsations if this pulsar had an X-ray luminosity and pulse shape similar to PSR B1937+21.

In long term timing of millisecond pulsars, for example with pulsar timing arrays like NANOGrav, it is very unusual for frequency derivatives above $\dot{f}$ to be significantly detected. One notable counterexample is PSR J1024$-$0719, which has a significantly detected second spin frequency derivative (Kaplan et al., 2016; Bassa et al., 2016b), ascribed to a long period orbit. A discussion of the measured timing noise and the possibility of PSR J1555$-$2908 being in a triple system with a long period outer orbit will be presented in Nieder et al. (2021, in prep.).

Another unusual trait of this MSP is the aligned radio and gamma-ray peaks, with similar morphology. These MSPs, dubbed Type A (for “aligned”) by Espinoza et al. (2013) or Class II by Johnson et al. (2014), account for only a few percent of the known gamma-ray MSPs in the Fermi Third Pulsar Catalog (3PC, in prep). This subset of MSPs are generally highly energetic, but the strongest correlation seems to be with the inferred magnetic field at the light cylinder ($B_{\mathrm{LC}}$; Espinoza et al. 2013). In the radio band they also generally have steeper spectra and lower levels of linear polarization than the general MSP population, and several are known to emit giant pulses. PSR J1555$-$2908 follows this trend, with its high $\dot{E}$, very large $B_{\mathrm{LC}}$, and steep ($\alpha=-2.5$; Frail et al. 2018) radio spectral index; however, we do not see evidence for giant pulses in the long 820 MHz GBT observation. Unfortunately, an instrumental issue with the GUPPI backend prevented us from making a measurement of the linear polarization fraction. The most successful models for these pulsars invoke co-located radio and gamma-ray emission regions in the outer magnetosphere, such as altitude-limited Two Pole Caustic models (alTPC; Johnson et al. 2014).

As has been the case for many of the MSPs discovered in searches of LAT sources, this pulsar is in an interacting binary system. The orbital solution gives a minimum companion mass of $0.052M_{\odot}$, for an assumed neutron star mass of 1.4 $M_{\odot}$. The very low mass companion [much lower than the white dwarf mass predicted by Tauris & Savonije (1999) of 0.17 $M_{\odot}$] and radio eclipses put it in the black widow class. As shown in Figure 6, the minimum companion mass is relatively large for MSPs in this orbital period range, which tend to have minimum companion masses of 0.02–0.03 $M_{\odot}$, but the eclipse around phase 0.25 indicates that the system is nearly edge-on, so the true companion mass is close to the minimum value. The eclipse spans $\approx 10\%$ of the orbit, which is twice the width of the companion’s Roche lobe (see below) in an edge-on system, consistent with absorption by an evaporated wind. The measured eccentricity is very low, but in family for black widow type systems.

Our optical observations confirm that the optical source detected by Gaia is the heated companion star to PSR J1555$-$2908. The amplitude of the modulation, likely more than 4 magnitudes from maximum to the non-detection at minimum, is among the largest seen in black widow companions (e.g., Breton et al., 2013; Draghis et al., 2019), which is consistent with heating being powered by the unusually high spin-down power. With a light curve in only one filter, and incomplete detection around the orbit, it is not possible to constrain a detailed model of the system geometry and irradiation of the companion star by the pulsar wind. Instead, we make a back-of-the-envelope calculation to check whether the peak observed flux of the companion is consistent with heating supplied by the spin-down power of the pulsar, and see if this constrains the temperature of the heated side of the companion. For an assumed isotropic pulsar wind and a Roche-lobe filling companion, the maximum irradiating power can be approximated as $L_{\rm irr}=\eta\,\dot{E}\,r_{L}^{2}/4a^{2}$, where $r_{L}$ is the radius of the Roche lobe of the companion, $a$ is the orbital separation, $\eta$ is the heating efficiency, and $\dot{E}\approx 3.1\times 10^{35}$ erg s^-1 for PSR J1555$-$2908. Assuming masses of $1.4\,M_{\odot}$ and $0.052\,M_{\odot}$, the ratio $r_{L}/a=0.154$ according to Eggleton (1983). The maximum irradiating power is then $L_{\rm irr}=1.8\times 10^{33}\eta$ erg s^-1.

In order to compare this to the observed optical flux at orbital phase 0.75, we assume an orbital inclination of $i\approx 90^{\circ}$, which is justified by the radio eclipses and the large amplitude of optical modulation around the orbit. We also correct the observed $r=20.4$ for $A_{r}=0.316$, the maximum extinction along the line of sight (Schlafly & Finkbeiner, 2011). Converting magnitude to flux density using $f^{\rm Vega}_{\lambda,\rm eff}=2.40\times 10^{-9}$ erg cm^-2 s^-1 Å^-1 (Fukugita et al., 1995), the absorption-corrected flux density at $\lambda^{\rm Vega}_{\rm eff}=6205$ Å is $f_{\lambda,{\rm eff}}=2.2\times 10^{-17}$ erg cm^-2 s^-1 Å^-1. For the DM distance we assume either 2.65 kpc or 7.55 kpc according to the two models of the electron density distribution referenced above.

At $d=2.65$ kpc, we find that the temperature of the heated half of the star is required to be $\geq 5,000$ K in order for the star to fit within its Roche lobe. At $T_{\rm h}=5,000$ K, the efficiency $\eta$ would be $0.05$. For higher temperature the heated area would decrease, allowing the star to partially fill its Roche lobe, and the efficiency would increase modestly. Assuming instead that $d=7.55$ kpc, the hot side must have $T_{\rm h}>9,000$ K to fit within the Roche-lobe, but the efficiency becomes $\eta\geq 0.44$. Any temperature above 14,000 K would require $\eta>1$. Thus, the range of plausible models would be much more limited at the larger distance.

We have additionally obtained 12 hr of simultaneous photometric observations in the $u$, $g$ and $i$ bands with the high-speed multi-band imager ULTRACAM on the 3.6m New Technology Telescope at ESO La Silla, as well as 0.65 orbits of optical spectroscopy, covering the optical maximum, with X-SHOOTER on ESO’s Very Large Telescope. These observations and their modelling will be presented in Kennedy et al. (2021, in prep.).