Discovery of a Gamma-ray Black Widow Pulsar by GPU-accelerated Einstein@Home
L. Nieder, C. J. Clark, D. Kandel, R. W. Romani, C. G. Bassa, B. Allen, A. Ashok, I. Cognard, H. Fehrmann, P. Freire, R. Karuppusamy, M. Kramer, D. Li, B. Machenschalk, Z. Pan, M. A. Papa, S. M. Ransom, P. S. Ray, J. Roy, P. Wang, J. Wu, C. Aulbert, E. D. Barr, B. Beheshtipour, O. Behnke, B. Bhattacharyya, R. P. Breton, F. Camilo, C. Choquet, V. S. Dhillon, E. C. Ferrara, L. Guillemot, J. W. T. Hessels, M. Kerr, S. A. Kwang, T. R. Marsh, M. B. Mickaliger, Z. Pleunis, H. J. Pletsch, M. S. E. Roberts, S. Sanpa-arsa, B. Steltner
DD RAFT VERSION O CTOBER
23, 2020Typeset using L A TEX twocolumn style in AASTeX63
Discovery of a Gamma-ray Black Widow Pulsar by GPU-accelerated Einstein@Home
L. N
IEDER ,
1, 2
C. J. C
LARK , D. K
ANDEL , R. W. R
OMANI , C. G. B
ASSA , B. A
LLEN ,
1, 6, 2
A. A
SHOK ,
1, 2
I. C
OGNARD ,
7, 8
H. F
EHRMANN ,
1, 2
P. F
REIRE , R. K
ARUPPUSAMY , M. K
RAMER ,
9, 3
D. L I ,
10, 11
B. M
ACHENSCHALK ,
1, 2
Z. P AN , M. A. P
APA ,
1, 6, 2
S. M. R
ANSOM , P. S. R AY , J. R OY , P. W
ANG , J. W U , C. A
ULBERT ,
1, 2
E. D. B
ARR , B. B
EHESHTIPOUR ,
1, 2
O. B
EHNKE ,
1, 2
B. B
HATTACHARYYA , R. P. B
RETON , F. C
AMILO , C. C
HOQUET , V. S. D
HILLON ,
17, 18
E. C. F
ERRARA ,
19, 20
L. G
UILLEMOT ,
7, 8
J. W. T. H
ESSELS ,
5, 21
M. K
ERR , S. A. K
WANG , T. R. M
ARSH , M. B. M
ICKALIGER , Z. P
LEUNIS ,
24, 25
H. J. P
LETSCH , M. S. E. R
OBERTS ,
26, 27
S. S
ANPA - ARSA , AND
B. S
TELTNER
1, 21
Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), 30167 Hannover, Germany Leibniz Universität Hannover, 30167 Hannover, Germany Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, M13 9PL, UK KIPAC/Dept. of Physics, Stanford University, Stanford, CA 94305, USA ASTRON, The Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands Department of Physics, University of Wisconsin–Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA Laboratoire de Physique et Chimie de l’Environnement et de l’Espace, Université d’Orléans / CNRS, F-45071 Orléans Cedex 02, France Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France Max-Planck-Institut für Radioastronomie, auf dem Hügel 69, 53121 Bonn, Germany National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China NAOC-UKZN Computational Astrophysics Centre, University of KwaZulu-Natal, Durban 4000, South Africa National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA USA 22903 Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA National Centre for Radio Astrophysics, Tata Institute of Fundamental Research, Pune 411 007, India South African Radio Astronomy Observatory, 2 Fir Street, Black River Park, Observatory 7925, South Africa Résidence Le Dauphiné, rue Jean Bleuzen, Vanves, France Department of Physics and Astronomy, University of Sheffield, Sheffield S3 7RH, UK Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA Department of Astronomy, University of Maryland, College Park, MD 20742, USA Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA Astronomy and Astrophysics Group, Department of Physics, University of Warwick, Coventry CV4 7AL, UK Department of Physics, McGill University, 3600 rue University, Montréal, QC H3A 2T8, Canada McGill Space Institute, McGill University, 3550 rue University, Montréal, QC H3A 2A7, Canada New York University Abu Dhabi, P.O. Box 129188, Abu Dhabi, UAE Eureka Scientific, Inc. 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA National Astronomical Research Institute of Thailand (Public Organization), 260 Moo 4, T. Donkaew, A. Maerim, Chiang Mai, 50180, Thailand (Received 2020 September 1; Revised 2020 September 22; Accepted 2020 September 25; Published 2020 October 22)
ABSTRACTWe report the discovery of 1 .
97 ms period gamma-ray pulsations from the 75 minute orbital-period bi-nary pulsar now named PSR J1653 − Fermi -Large Area Telescope gamma-ray source4FGL J1653.6 − Einstein@Home distributed volunteer computing sys-
Corresponding author: L. [email protected] a r X i v : . [ a s t r o - ph . H E ] O c t N IEDER ET AL .tem. The multi-dimensional parameter space was bounded by positional and orbital constraints obtained fromthe optical counterpart. More sensitive analyses of archival and new radio data using knowledge of the pulsartiming solution yield very stringent upper limits on radio emission. Any radio emission is thus either exception-ally weak, or eclipsed for a large fraction of the time. The pulsar has one of the three lowest inferred surfacemagnetic-field strengths of any known pulsar with B surf ≈ × G. The resulting mass function, combined withmodels of the companion star’s optical light curve and spectra, suggests a pulsar mass (cid:38) M (cid:12) . The companionis light-weight with mass ∼ . M (cid:12) , and the orbital period is the shortest known for any rotation-poweredbinary pulsar. This discovery demonstrates the Fermi -Large Area Telescope’s potential to discover extremepulsars that would otherwise remain undetected.
Keywords: gamma rays: stars — pulsars: individual (PSR J1653 − INTRODUCTIONThe
Fermi
Large Area Telescope (LAT) source4FGL J1653.6 − − − − Einstein@Home . Such searches are very compu-tationally demanding, and would take decades to centurieson a single computer while still taking weeks or months on
Einstein@Home . Thus, the search methods are specificallydesigned to ensure efficiency (Nieder et al. 2020). One keyelement is the use of constraints derived from optical obser-vations. The companion’s pulsar-facing side is heated by thepulsar wind, leading to a periodically varying optical lightcurve. This permits the orbital period P orb and other orbitalparameters to be tightly constrained (for a feasible search theuncertainty ∆ P orb needs to be less than a few milliseconds).In addition, because the sky position of the optical source istypically known to high precision (sub-milliarcsecond level),a search over position parameters is not needed.Here we present the discovery and analysis of gamma-raypulsations from PSR J1653 − − GAMMA-RAY PULSATIONS2.1.
Data preparation
We searched for gamma-ray pulsations in the arrival timesof photons observed by the
Fermi
LAT (Atwood et al. 2009)between 2008 August 3 and 2018 April 16 (MJDs 54 , , SOURCE -class photons accordingto the
P8R2_SOURCE_V6 (Atwood et al. 2013) instrumentresponse functions (IRFs) , with reconstructed incidence an-gles within a 5 ◦ region of interest (RoI) around the puta-tive pulsar position, energies above 100 MeV, and zenith an-gles below 90 ◦ . Here, we used the presumptive compan-ion’s position as reported in the Gaia
DR2 Catalog (here-after
Gaia catalog; Gaia Collaboration et al. 2018). The ce-lestial parameters (J2000.0) are α = 16 h m . s δ = − ◦ (cid:48) . (cid:48)(cid:48) σ uncertainties on the last dig-its reported in parentheses.Using the photon incidence angles and energies, we con-structed a probability or weight for each photon, w j ∈ [0 , j labels the photon: w j is the probability that the j thphoton originated from the posited source, as opposed to afore- or background source. These weights were computedby gtsrcprob , using the preliminary Fermi -LAT 8 yearsource catalog as a model for the flux within the RoI with-out performing a full spectral fit. Weighting the contribu-tion of each photon to a detection statistic in this way greatlyincreases the search sensitivity (Kerr 2011), and the distri-bution of weights can be used to predict expected signal-to-noise ratios (Nieder et al. 2020). See https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_essentials.html https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/ ISCOVERY OF
PSR J1653 − THROUGH G AMMA R AYS N = 354 ,
009 photons,collected over a period of 3 ,
542 days. The properties of thedetection statistics (semicoherent power S , coherent power P , and H statistic) depend upon the lowest moments of theweights, which are N (cid:88) j =1 w j ≈ , N (cid:88) j =1 w j ≈ , and N (cid:88) j =1 w j ≈ . These moments determine the ultimate sensitivity to a partic-ular pulse profile and pulsed fraction, as given in Eq. (11) inNieder et al. (2020).Following the pulsar discovery, we extended this datasetto 2020 February 23 (MJD 58 , P8R3_SOURCE_V2
IRFs (Bruel et al. 2018), a larger maxi-mum zenith angle of 105 ◦ , and using the Fermi -LAT FourthSource Catalog (hereafter 4FGL; Abdollahi et al. 2020) asthe RoI model for the photon probability weight computa-tions. 2.2.
Search
The binary-pulsar search methods are described by Niederet al. (2020), which are a generalization and extension of theisolated-pulsar search methods from Pletsch & Clark (2014).The searched ranges are guided by the known millisecondpulsar (MSP) population in the Australia Telescope NationalFacility (ATNF) Pulsar Catalogue (Manchester et al. 2005).For the spin frequency, we searched f ∈ [0 , . Thespin-frequency derivative was expected to be in the range ˙ f ∈ [ − − ,
0] Hz s − .The sky position of the candidate optical counterpart isconstrained to high precision in the Gaia catalog, so no as-trometric search is required. The proper motion measured by
Gaia for the optical counterpart was ignored for the search.2.2.1.
Orbital Constraints from Optical Observations
The orbital-period estimate of Romani et al. (2014) wasderived from Southern Astrophysical Research (SOAR),WIYN, and Catalina Sky Survey (CSS) observations. Thesewere augmented by new 350s SOAR Goodman HighThroughput Spectrograph (GHTS) g (cid:48) , r (cid:48) , i (cid:48) exposures (63 g (cid:48) ,75 r (cid:48) , 42 i (cid:48) ) from MJD 56 , .
074 – 56 , . g (cid:48) , r (cid:48) , and i (cid:48) exposures obtained by Kong et al.(2014) using the Wide Field camera (WFC) on the 2.5m IsaacNewton Telescope (INT) on La Palma. For these two datasets, the scatter about the light-curve trends was apprecia-bly larger than the very small statistical errors; we thus add The upper limit has been chosen to be sensitive to pulsars spinning at up to750 Hz, which have two-peaked pulse profiles where the peaks are half arotation apart (see also Pletsch & Clark 2014). Note that the current recordspin frequency is 716 Hz (Hessels et al. 2006). .
03 mag in quadrature to account for unmodeled fast vari-ability and/or photometry systematics. To further refine theorbital-period uncertainty, we obtained additional observa-tions in u (cid:48) , g (cid:48) , and i (cid:48) using the high-speed multi-band im-ager ULTRACAM (Dhillon et al. 2007) on the 4.2m WilliamHerschel Telescope (WHT) on two nights (MJDs 57 ,
170 and57 , > g (cid:48) , 151 r (cid:48) , 45 i (cid:48) ) onMJD 57 ,
988 – 57 , g (cid:48) , r (cid:48) , i (cid:48) filter fluxes were refer-enced to in-field PanSTARRS catalog sources, and then con-verted to the Sloan Digital Sky Survey (SDSS) scale. The u (cid:48) photometry was calibrated against an SDSS standard star ob-served on MJD 57 , ∼ .
05 mag systematicuncertainties in g (cid:48) , r (cid:48) , and i (cid:48) , with uncertainties as large as ∼ . u (cid:48) .We constrained the orbital period using the multi-bandLomb-Scargle periodogram method (VanderPlas & Ivezi´c2015, excluding the u (cid:48) ULTRACAM data, as the modulationhas very low signal-to-noise ratio in this band). To infer rea-sonable statistical uncertainties, we fit for and removed con-stant magnitude offsets, consistent with our estimated cali-bration uncertainties, between each night’s observations ineach band, and additionally rescaled the magnitude uncer-tainties to obtain a reduced chi-square of unity. This con-strained the orbital period to P orb = 0 . ± . × − days, where the quoted uncertainty is the 1 σ statisticaluncertainty. For the pulsation search, we chose to search the3 σ range around this value.In Romani et al. (2014), the time of the pulsar’s ascend-ing node, T asc , was estimated from the photometric lightcurve. However, the optical maximum is distinctly asym-metric (see Section 3.1), which can bias orbital phase esti-mates. We therefore used the spectroscopic radial-velocitymeasurements from Romani et al. (2014), folded at the or-bital period obtained above, and fit the phase of a sinusoidalradial-velocity curve, finding T asc = MJD 56 , . ± . × − . However, as radial velocities may still be slightlybiased by asymmetric heating, we elected to search a widerange around this value, corresponding to ± σ .For the projected semimajor-axis parameter x = a sin i / c ,we decided to start searching x ∈ [0 , .
1] s, with the intentionto go to larger values in the case of no detection. For a pulsarmass of 1 . M (cid:12) , this would cover the companion mass rangeup to 0 . M (cid:12) and would include companion masses of allknown “black-widow” systems as well as some of the lower-mass “redback” systems (Roberts 2013; Strader et al. 2019). N IEDER ET AL .Here, a is the pulsar’s semimajor axis, i denotes the inclina-tion angle, and c is the speed of light. As described in Niederet al. (2020), we expected x ∈ [0 , .
2] s based on the com-panion’s velocity amplitude reported by Romani et al. (2014)and the masses expected for “spider” companions, i.e. black-widow or redback companions.2.2.2.
Search grids
To cover the relevant orbital-parameter space in { x , P orb , T asc } , we use optimized grids (Fehrmann & Pletsch2014). These grids use as few points as possible still ensur-ing that a signal within the relevant space should be detected.Furthermore, they are able to cover the orbital-parameterspace efficiently even though the required density dependson one of the orbital parameters, x .Key to building an optimized grid is to know how thesignal-to-noise ratio drops due to offsets from the true pul-sar parameters. This is estimated using a distance metric onthe orbital-parameter space (Nieder et al. 2020). In our case,the three-dimensional grid was designed to have a worst-casemismatch ¯ m = 0 .
2, i.e. not more than 20% of the (semicoher-ent or coherent) signal power should be lost due to orbital-parameter offsets. Of most relevance is that 99% of ran-domly injected orbital-parameter points have a mismatch be-low ¯ m = 0 .
04 to the closest grid point.Due to the f -dependency of the required grid-point den-sity, we search f in steps, and build the corresponding orbitalgrids prior to the start of the search on the computing clusterATLAS in Hannover (Aulbert & Fehrmann 2008).2.2.3. Einstein@Home
Searching the 5-dimensional parameter space { f , ˙ f , x , P orb , T asc } is a huge computational task with over10 trials. Thus, the first (computing-intensive) searchstages were performed on Einstein@Home , a distributed vol-unteer computing system (Allen et al. 2013). As done forradio pulsar searches previously, the search code utilizes theapproximately 10 ,
000 GPUs active on
Einstein@Home fora computing speed-up of ∼
10, comparing the runtimes onCPUs and GPUs.The parameter space is divided into more than one mil-lion regions. Searching one of these is called a “work unit”.These work units are sent to computers participating in
Ein-stein@Home , and are searched when the computer is other-wise idle. Depending on the system, searching a work unittakes between half an hour and up to a few hours of compu-tational time. In total, the search would have taken more than50 years on a single computer, but using
Einstein@Home ittook less than two weeks.2.2.4.
Gamma-ray detection
The search process involves multiple stages in which semi-coherent statistics are constructed, and the most significant candidates are passed on to fully coherent follow-up stages(for full details of the search pipeline and signal-to-noise ra-tio definitions, see Nieder et al. 2020). In the last semicoher-ent stage, a candidate found at a frequency of 1016 Hz hadsignal-to-noise ratio S = 8 .
6, which we now associate withPSR J1653 − P / P ≈ .
97 ms (or f ≈
508 Hz), while the actualsearch revealed an alias at twice the pulsar frequency. Thismay be because the signal has significant power in the secondharmonic.Note that the signal was found outside the 3 σ range in T asc from the constraints reported in this work, and outside the 3 σ range given by Romani et al. (2014). This can be caused byasymmetric heating (see Section 2.2.1).2.3. Timing
The parameters used in the phase model to describe thepulsar’s rotation are measured in a timing analysis. We usethe timing methods as explained in Clark et al. (2017), whichare an extension of the methods by Kerr et al. (2015). Thebasic principle is that the parameter space around the dis-covery parameters is explored using a Monte Carlo samplingalgorithm with a template pulse profile.To marginalize over the pulse-profile template, we vary thetemplate parameters as described in Nieder et al. (2019). Inthe case of PSR J1653 − ˙ P is one of the lowest of allknown pulsars. To estimate the intrinsic ˙ P we account for theShklovskii effect (Shklovskii 1970), and the Galactic acceler-ation (see, e.g., Damour & Taylor 1991). The results are sum-marized in Table 1. The observed contribution due to the dif-ference in Galactic acceleration of the Sun and the pulsar iscomputed with R Sun = 8 .
21 kpc, z Sun = 14 pc, and the Galacticpotential model
PJM17_best.Tpot (McMillan 2017), asimplemented in their code . For PSR J1653 − R J1653 = 7 .
48 kpc, and z J1653 = 367 pc, assuming d = 840 pc(see Table 2). The contributions parallel and perpendicular tothe Galactic disk nearly cancel each other, so that the choice https://github.com/PaulMcMillan-Astro/GalPot ISCOVERY OF
PSR J1653 − THROUGH G AMMA R AYS Table 1.
Timing solution for PSR J1653 − . Gaia catalogR.A., α (J2000.0) 16 h m . s δ (J2000.0) − ◦ (cid:48) . (cid:48)(cid:48) . µ α cos δ (mas yr − ) − . ± . µ δ (mas yr − ) − . ± . a , ϖ (mas) 1 . ± . f (Hz) 508 . ˙ f (Hz s − ) − . × − Spin period, P (ms) 1 . ˙ P (s s − ) 2 . × − Proj. semimajor axis, x (s) 0 . P orb (days) 0 . T asc (MJD) 56513 . d = 840 pcShklovskii spin-down, ˙ P Shk (s s − ) 1 . × − Galactic acceleration spin-down, ˙ P Gal (s s − ) − . × − Spin-down power, ˙ E (erg s − ) 4 . × Surface B -field, B surf (G) 4 . × Light-cylinder B -field, B LC (G) 5 . × Characteristic age, τ c (Gyr) 37Gamma-ray luminosity b , L γ (erg s − ) 2 . × Gamma-ray efficiency, n γ = L γ / ˙ E . OTE —The JPL DE405 solar system ephemeris has been used, andtimes refer to TDB. a Corresponds to a model-independent distance d = 533 + − pc, but forthe derived parameters the consistent distance d = 840 + − pc derivedfrom optical modeling is used (see Table 2). b Taken from 4FGL Source Catalog (Abdollahi et al. 2020). of the potential and its relevant parameters have a seeminglylarge effect on the actual small value of ˙ P Gal , and can evenchange the sign. However, the overall kinematic contributionto the observed ˙ P is dominated by the Shklovskii term, andits uncertainty by the uncertainty in the distance estimate.The estimated intrinsic spin-down is ˙ P int = 8 . × − s s − for distance d = 840 pc. Figure 1.
Integrated pulse-profile and phase-time diagram ofPSR J1653 − − ,
600 and 57 ,
000 when theLAT pointed more often toward the Galactic center.3.
MULTIWAVELENGTH AND MULTIMESSENGER3.1.
Optical Light-curve Modeling and System Masses
By modeling the optical light curves and radial velocitieswe can constrain the binary mass and distance and the sys-tem viewing angle. Comparing the individual filters betweennights suggest small δ m ≈ .
05 shifts in zero-points, consis-tent with the systematic estimates above. Correcting to matchthe individual filters, we then re-binned the light curve, plac-ing the photometry on a regular grid with points spaced by δφ = 0 . Lightkurve ; af-ter excision of a few obviously discrepant points, we retain N
IEDER ET AL . .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . Binary Phase . . . . . . . . F l u x ( e r g s − H z − c m − ) × − Model i r g u Figure 2. u (cid:48) , g (cid:48) , r (cid:48) , and i (cid:48) light curves for PSR J1653 − u (cid:48) , 239 g (cid:48) , 220 r (cid:48) and 245 i (cid:48) points for light-curve fitting(Fig. 2). This fitting is done with a version of the Icarus code of Breton et al. (2013) modified to include the effectof hot spots on the companion surface, likely generated byprecipitation of particles from the intrabinary shock (IBS) tocompanion magnetic poles (Sanchez & Romani 2017). Allparameter values and errors are determined by Markov ChainMonte Carlo (MCMC) modeling.The very shallow modulation of these light curves mightnormally be interpreted as indicating a small inclination i . However given the large companion radial-velocity am-plitude K = 666 . ± . − , implying a mass function f ( M ) = 1 . ± . M (cid:12) , measured by Romani et al. (2014),a small inclination would give an unphysical, large neutronstar mass. As noted in that paper, the light curves and spectrashow that a strong blue non-thermal veiling flux dominates atorbital minimum. With increasingly shallow modulation forthe bluer colors, this is also evident in the present photom-etry. Thus, the minimal model for this pulsar must includea non-thermal veiling flux. Although this is likely associ-ated with the IBS, we model it here as a simple power lawwith form f ν = f A ( ν/ Hz) − p . This flux is nearly constantthrough the orbit, although there are hints of phase structure,e.g. in r (cid:48) and i (cid:48) at φ B = 0 .
72 (see Fig. 2). Any model with-out such a power-law component is completely unacceptable.These fits prefer an A V slightly higher than, but consistentwith, the maximum in this direction (obtained by ∼
300 pc;Green et al. 2019) . https://doi.org/10.7910/DVN/2EJ9TX Table 2.
Light-curve fit results for PSR J1653 − i (deg) 79 . + . − . . + . − . Filling factor, f c . + . − . . + . − . Heating luminosity, L P (10 erg s − ) 3 . + . − . . + . − . Night-side temperature, T N (K) 3250 + − + − V -band extinction, A V . + . − . . + . − . Distance, d (pc) 830 + − + − Veiling flux norm, f A ( µ Jy) 101 . + . − . . + . − . Veiling flux index, p . + . − . . + . − . Spot azimuth, θ c (deg) ... 286 . + . − . Spot co-latitude, φ c (deg) ... − . + . − . Gaussian spot width, σ c (deg) ... 25 . + . − . Spot temperature increase, A c ... 0 . + . − . Neutron star mass, M NS ( M (cid:12) ) 1 . + . − . . + . − . Companion mass, M c ( M (cid:12) ) 0 . + . − . . + . − . χ / DoF 1 .
72 1 . OTE —Parameters from the best-fit light-curve/radial-velocitymodels, with and without a surface hot spot, including MCMCerrors.
In Fig. 2, one notices that the orbital maximum is slightlydelayed from φ B = 0 .
75, especially in the bluer colors. Suchasymmetric heating is most easily modeled adding a polarhot spot with location ( θ c , φ c ) and local temperature increase A c in a Gaussian pattern of width σ c ; when we include such acomponent, the fit improves greatly, with ∆ χ / DoF = − . − level, despite the extra degrees of freedom.We give the fit parameters for both models in Table 2. Notethat with the fine structure near maximum, the model is notyet fully acceptable ( χ / DoF ∼ . f γ =3 . × − erg cm − s − between 100 MeV and 100 GeV,our fit distance gives an isotropic gamma-ray luminosity L γ = 3 × erg s − , in good agreement with the L γ ≈ (10 erg s − ˙ E ) / heuristic luminosity law (Abdo et al. 2013), ISCOVERY OF
PSR J1653 − THROUGH G AMMA R AYS ˙ E . This luminosity isconsistent with the model for direct radiative heating of thecompanion. (2) Our fit distance is also consistent with themodel-independent, but lower-accuracy, distance from the Gaia parallax. Thus, the 840 pc distance seems reliable, al-though systematic effects probably dominate over the rathersmall ∼
50 pc statistical errors.Armed with the fits, we can estimate the companionmasses, correcting the observed radial-velocity amplitude(fit with a K-star template) for the temperature-dependentweighting of the absorption lines across the companion faceas in Kandel & Romani (2020). The results indicate substan-tial mass accretion, as expected for these ultra-short-periodsystems. With the preferred Veiled+HS model the mass sig-nificantly exceeds 2 . M (cid:12) , adding to the growing list of spi-der binaries in this mass range. Note that the inclination i uncertainty dominates the error in this mass determination.Broader range photometric studies, with better constraint onthe heating pattern, can reduce the i uncertainty.3.2. Radio pulsation searches
The pulsar position has been observed in radio multipletimes. Several searches were performed before the gamma-ray pulsation discovery, and a few very sensitive follow-upsearches afterward. Despite the more than 20 observationswith eight of the most sensitive radio telescopes, no radiopulsations have been found.The results of the radio searches are given in Table 3. Ob-servations are spread over 11 years, with observing frequen-cies ranging from 100 MHz up to 5 GHz. All orbital phaseshave been covered by most of the telescopes. Since therewas no detection, the table also gives upper limits derivedfrom the observations. For all but LOFAR, the data (botharchival and recent) were folded with the gamma-ray-derivedephemeris, and searched only over dispersion measure.The strictest upper limits on pulsed radio emission are8 µ Jy at 1 . µ Jy at 4 . µ Jy that Abdo et al. (2013) use todefine a pulsar to be “radio-quiet”. Note, that for the calcu-lation of the limits we included the parts of the orbit whereeclipses might be expected for spider pulsars. Thus, the limitconstrains the maximum emission of the system, and not themaximum emission from the pulsar alone.3.3.
Continuous gravitational waves
We search for nearly monochromatic, continuousgravitational waves (GWs) from PSR J1653 − and second observing runs of the Ad-vanced LIGO detectors (The LIGO Scientific Collaboration https://doi.org/10.7935/K57P8W9D https://doi.org/10.7935/CA75-FM95 et al. 2019). We assume that GWs are emitted at the first andsecond harmonic of the neutron star’s rotational frequency, aswould occur if the spin axis is misaligned with the principalaxes of the moment of inertia tensor (Jones 2010, 2015).We employ two different analysis procedures, which yieldconsistent results. The first is frequentist, based on the multi-detector maximum-likelihood F -statistic introduced by Cut-ler & Schutz (2005). The second is the Bayesian time-domain method (Dupuis & Woan 2005) as detailed by Pitkinet al. (2017), with triaxial non-aligned priors (Pitkin et al.2015). Both methods coherently combine data from the twodetectors, taking into account their antenna patterns and theGW polarization. The F -statistic search excludes data takenduring times when the relevant frequency bands are exces-sively noisy.The results are consistent with no GW emission. At twicethe rotation frequency, the F -statistic 95% confidence upperlimit on the intrinsic GW amplitude h is 4 . × − . The95% credible interval upper limit from the Bayesian analysison h = 2 C is 3 . × − . At the rotation frequency (onlychecked with the Bayesian method) the 95% confidence up-per limit on the amplitude C is 6 . × − .Since the dominant GW frequency might be mismatchedfrom twice the rotation frequency (Abbott et al. 2019a), weperformed an F -statistic search in a ± ˙ f -range. This yields larger upper limits on h , with mean value of 1 . × − in 10 mHz-wide bands.Full details are given in the supplementary materials.Our upper limits on h at twice the rotation frequencymay also be expressed as upper limits on the ellipticity (cid:15) ofthe pulsar (Abbott et al. 2019b). This is (cid:15) = 3 . × − × ( h / × − ) × (10 g cm / I zz ) × (840 pc / d ), where I zz isthe moment of inertia about the spin axis, and d is the dis-tance.As is the case for most known pulsars, it is unlikely that oursearches would have detected a GW signal. In fact, supposethat all of the rotational kinetic-energy losses associated withthe intrinsic spin-down are via GW emission. Then assum-ing the canonical I zz = 10 g cm , this would imply a “spin-down” ellipticity (cid:15) sd = 4 . × − , which is a factor ∼ DISCUSSION AND CONCLUSIONSPSR J1653 − f = 508 Hz. The 75 min orbital periodis shorter than for any other known rotation-powered pul-sar, with the previous record being PSR J1311 − IEDER ET AL . Table 3.
Summary of radio searches for PSR J1653 − µ Jy) Reference / SurveyEffelsberg 1210–1510 2010 May 26, 21:33 1920 0 . .
31 63 Barr et al. (2013)Effelsberg 1210–1510 2014 Aug 26, 20:27 4600 0 . .
17 41Effelsberg 4608–5108 2014 Aug 29, 18:52 4600 0 . .
65 33Effelsberg 4608–5108 2020 Jun 18, 22:09 11820 0 . .
48 20FAST 1050–1450 2020 Jun 04, 16:30 2036 0 . .
25 8 Li et al. (2018)GBT 720–920 2009 Sep 20, 00:49 3200 0 . .
65 51GBT 720–920 2010 Dec 13, 21:04 1300 0 . .
20 80GBT 720–920 2011 Dec 22, 12:11 2400 0 . .
27 59 Sanpa-arsa (2016)GBT 305–395 2012 Feb 22, 14:31 1700 0 . .
65 301GBT 1700–2300 2014 Nov 18, 14:28 1200 0 . .
63 43GBT 1700–2300 2014 Nov 20, 13:56 2400 0 . .
98 30GBT 1700–2300 2014 Nov 21, 22:38 1800 0 . .
07 35GBT 720–920 2017 Jan 28, 13:20 1200 0 . .
24 83GMRT 591–623 2011 Feb 02, 02:32 1800 0 . .
34 730 Bhattacharyya et al.GMRT 306–338 2012 May 15, 22:31 1800 0 . .
06 990 (2013, 2020, in prep.)GMRT 306–338 2012 Jun 11, 17:49 1800 0 . .
95 990 "GMRT 591–623 2014 Aug 19, 13:44 1800 0 . .
54 270 "GMRT 591–623 2014 Aug 30, 11:17 1800 0 . .
38 270 "GMRT 591–623 2015 Dec 28, 03:55 1800 0 . .
13 270 "LOFAR 110–180 2017 Mar 15, 04:18 15 ×
320 Full orbit 6 ,
200 Bassa et al. (2017)LOFAR 110–180 2017 Apr 15, 02:20 15 ×
320 Full orbit 6 ,
200 "Lovell 1332–1732 2019 Mar 15, 01:34 5400 0 . .
77 82Lovell 1332–1732 2019 Mar 16, 02:53 5400 0 . .
08 82Lovell 1332–1732 2019 Mar 17, 01:47 5400 0 . .
45 82Nançay 1230–1742 2014 Aug 20, 18:33 1850 0 . .
53 77 Desvignes et al. (2013)Parkes 1241–1497 2016 Nov 05, 06:17 3586 0 . .
06 178 Camilo et al. (2016)N
OTE —The columns show the telescope used, the observed frequency range, the start time and data span, the range of orbitalphases covered, the resulting limit on a pulsed component, and a reference with relevant details. The orbital phase is given inorbits, and ranges > − at the lowest frequencies to DM = 350 pc cm − at the highest frequencies.To estimate the limit on the pulsed component, we used Eq. (6) from Ray et al. (2011) assuming a pulse width of 0 . P , and athreshold signal-to-noise ratio S/N min = 7. a 93 min orbit (Pletsch et al. 2012). The inferred surfacemagnetic field is possibly the weakest, depending on theShklovskii correction.The discovery was enabled by constraints on the sky-position and orbital parameters from optical observations,together with efficient search techniques and the large com-puting power of the distributed volunteer computing system Einstein@Home . The detection proves that the optically vari-able candidate counterpart (Kong et al. 2014; Romani et al.2014) is indeed the black-widow-type binary companion toPSR J1653 − − Gaia measurements of the parallax, ϖ = 1 . ± .
01 mas, imply a distance d = 530 + − pc. A con-sistent, but tighter constraint is given by our optical modelingwith d = 840 + − pc. The proper motion (see Table 1) is alsomeasured with good precision ( Gaia and our timing are inagreement).PSR J1653 − ˙ P = 2 . × − s s − ).The intrinsic ˙ P = 8 . × − s s − (accounting for Galactic ac-celeration and Shklovskii effects) is even smaller. In Fig. 3,PSR J1653 − P - ˙ P diagram, alongside theknown radio and gamma-ray pulsar population outside ofglobular clusters. ISCOVERY OF
PSR J1653 − THROUGH G AMMA R AYS − − − Spin period P [s] − − − − − − − − Sp i n - p e r i o dd e r i v a t i v e ˙ P [ ss − ] k y r M y r G y r e r g s − e r g s − e r g s − G G G G Figure 3.
Newly detected PSR J1653 − P – ˙ P diagramof the known pulsar population outside of globular clusters. TheMSP population is shown magnified in the inset. LAT pulsars aremarked in green (isolated by a cross and binary by a circle). Non-LAT pulsars in the ATNF are marked in gray (isolated by a plusand binary by a square). The lines show constant surface magnetic-field strength (dashed-dotted), characteristic age (dotted), and spin-down power (dashed). The spin period and intrinsic spin-periodderivative of PSR J1653 − The intrinsic ˙ P can be used to estimate the pulsar’sspin-down power ˙ E , surface magnetic-field strength B surf ,magnetic-field strength at the light cylinder B LC , and char-acteristic age τ c . These are given in Table 1 for d = 840 pc.Constant lines of ˙ E , B surf , and τ c are displayed in Fig. 3 toshow the distance-dependent ranges.Spider pulsars in very-short-period orbits are difficult todiscover with traditional radio searches. Even though wecan now fold the radio data with the exact parameters,PSR J1653 − − − − is very high, assuming a filled Roche lobe (Eggleton 1983).Using the filling factor from optical modeling, the averagecompanion density 73 g cm − is even higher. The high den-sity and the compact orbit suggest that the companion may bea helium white-dwarf remnant, and that the system may haveevolved from an ultracompact X-ray binary (Sengar et al.2017; Kaplan et al. 2018). In addition, simulations predictevolved ultracompact X-ray binaries to have orbital periodsof around 70 −
80 min (van Haaften et al. 2012), consistentwith the 75 min orbital period from PSR J1653 − − − IEDER ET AL .operated on the island of La Palma by the Isaac NewtonGroup of Telescopes in the Spanish Observatorio del Roquede los Muchachos of the Instituto de Astrofísica de Canarias.Based on observations made with the Isaac Newton Tele-scope (program I17BN005) operated on the island of LaPalma by the Isaac Newton Group of Telescopes in the Span-ish Observatorio del Roque de los Muchachos of the Institutode Astrofísica de Canarias. This paper makes use of data ob-tained from the Isaac Newton Group of Telescopes Archivewhich is maintained as part of the CASU Astronomical DataCentre at the Institute of Astronomy, Cambridge.We acknowledge support of the Department of Atomic En-ergy, Government of India, under project No. 12-R&D-TFR-5.02-0700 for the GMRT observations. The GMRT is run bythe National Centre for Radio Astrophysics of the Tata In-stitute of Fundamental Research, India. The Nançay RadioObservatory is operated by the Paris Observatory, associatedwith the French Centre National de la Recherche Scientifique(CNRS). We acknowledge financial support from the “Pro-gramme National Hautes Energies” (PNHE) of CNRS/INSU,France. This Letter is based (in part) on data obtained withthe International LOFAR Telescope (ILT) under project codeLC7_018. LOFAR (van Haarlem et al. 2013) is the LowFrequency Array designed and constructed by ASTRON.The National Radio Astronomy Observatory is a facility ofthe National Science Foundation operated under cooperativeagreement by Associated Universities, Inc. The Green BankObservatory is a facility of the National Science Foundationoperated under cooperative agreement by Associated Univer-sities, Inc. FAST is a Chinese national mega-science facility,built and operated by NAOC. Partly based on observationswith the 100 m telescope of the MPIfR (Max-Planck-Institutfür Radioastronomie) at Effelsberg.The
Fermi -LAT Collaboration acknowledges generous on-going support from a number of agencies and institutes thathave supported both the development and the operation ofthe LAT as well as scientific data analysis. These include theNational Aeronautics and Space Administration and the De-partment of Energy in the United States, the Commissariat àl’Energie Atomique and the Centre National de la RechercheScientifique/Institut National de Physique Nucléaire et dePhysique des Particules in France, the Agenzia Spaziale Ital- iana and the Istituto Nazionale di Fisica Nucleare in Italy, theMinistry of Education, Culture, Sports, Science and Tech-nology (MEXT), High Energy Accelerator Research Orga-nization (KEK) and Japan Aerospace Exploration Agency(JAXA) in Japan, and the K. A. Wallenberg Foundation, theSwedish Research Council, and the Swedish National SpaceBoard in Sweden.Additional support for science analysis during the oper-ations phase is gratefully acknowledged from the IstitutoNazionale di Astrofisica in Italy and the Centre Nationald’Études Spatiales in France. This work performed in partunder DOE Contract DE-AC02-76SF00515.The authors thank the LIGO Scientific Collaboration foraccess to the data and gratefully acknowledge the support ofthe United States National Science Foundation (NSF) for theconstruction and operation of the LIGO Laboratory and Ad-vanced LIGO as well as the Science and Technology Facil-ities Council (STFC) of the United Kingdom, and the Max-Planck-Society (MPS) for support of the construction of Ad-vanced LIGO. Additional support for Advanced LIGO wasprovided by the Australian Research Council. This researchhas made use of data, software, and/or web tools obtainedfrom the LIGO Open Science Center (https://losc.ligo.org),a service of LIGO Laboratory, the LIGO Scientific Collab-oration and the Virgo Collaboration, to which the authorshave also contributed. LIGO is funded by the U.S. NationalScience Foundation. Virgo is funded by the French CentreNational de Recherche Scientifique (CNRS), the Italian Isti-tuto Nazionale della Fisica Nucleare (INFN), and the DutchNikhef, with contributions by Polish and Hungarian insti-tutes.
Software:
Fermi
Science Tools,
MultiNest (Ferozet al. 2019), ULTRACAM software pipelines,
Icarus (Breton et al. 2012), psrqpy (Manchester et al. 2005;Pitkin 2018),
Astropy (Astropy Collaboration et al. 2013,2018), matplotlib (Hunter 2007),
NumPy (Oliphant2006; van der Walt et al. 2011),
GalPot (McMillan 2017),
Lightkurve (Lightkurve Collaboration et al. 2018),
PRESTO (Ransom et al. 2002), LALSuite (LIGO ScientificCollaboration 2018)REFERENCES
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, PhRvD, 99,122002, doi: 10.1103/PhysRevD.99.122002—. 2019b, ApJ, 879, 10, doi: 10.3847/1538-4357/ab20cbAbdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJS, 183,46, doi: 10.1088/0067-0049/183/1/46Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17,doi: 10.1088/0067-0049/208/2/17 Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247,33, doi: 10.3847/1538-4365/ab6bcbAllen, B., Knispel, B., Cordes, J. M., et al. 2013, ApJ, 773, 91,doi: 10.1088/0004-637X/773/2/91Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science,324, 1411, doi: 10.1126/science.1172740
ISCOVERY OF
PSR J1653 − THROUGH G AMMA R AYS IEDER ET AL ..