Discovery of a Gamma-Ray Black Widow Pulsar by GPU-accelerated Einstein@Home

We report the discovery of 1.97 ms period gamma-ray pulsations from the 75 minute orbital-period binary pulsar now named PSR J1653−0158. The associated Fermi Large Area Telescope gamma-ray source 4FGL J1653.6−0158 has long been expected to harbor a binary millisecond pulsar. Despite the pulsar-like gamma-ray spectrum and candidate optical/X-ray associations—whose periodic brightness modulations suggested an orbit—no radio pulsations had been found in many searches. The pulsar was discovered by directly searching the gamma-ray data using the GPU-accelerated Einstein@Home distributed volunteer computing system. The multidimensional parameter space was bounded by positional and orbital constraints obtained from the optical counterpart. More sensitive analyses of archival and new radio data using knowledge of the pulsar timing solution yield very stringent upper limits on radio emission. Any radio emission is thus either exceptionally weak, or eclipsed for a large fraction of the time. The pulsar has one of the three lowest inferred surface magnetic-field strengths of any known pulsar with Bsurf ≈ 4 × 107 G. The resulting mass function, combined with models of the companion star’s optical light curve and spectra, suggests a pulsar mass ≳2 M⊙. The companion is lightweight with mass ∼0.01 M⊙, and the orbital period is the shortest known for any rotation-powered binary pulsar. This discovery demonstrates the Fermi Large Area Telescope's potential to discover extreme pulsars that would otherwise remain undetected.


INTRODUCTION
The Fermi Large Area Telescope (LAT) source 4FGL J1653.6−0158 is a bright gamma-ray source, and the brightest remaining unassociated source (Saz Parkinson et al. 2016).It was first seen by the Energetic Gamma Ray Experiment Telescope (EGRET, Hartman et al. 1999), and was also listed in the LAT Bright Gamma-ray source list (Abdo et al. 2009) more than a decade ago.While pulsars were discovered in several other sources from this list (see, e.g., Ransom et al. 2011), the origin of 4FGL J1653.6−0158remained unidentified.The detection of a variable X-ray and optical candidate counterpart with 75 min-period consistent with the gamma-ray position of 4FGL J1653.6−0158provided strong evidence of it being a binary gamma-ray pulsar (Kong et al. 2014;Romani et al. 2014).
To identify the neutron star in 4FGL J1653.6−0158,we carried out a binary-pulsar search of the gamma rays, using the powerful GPU-accelerated distributed volunteer computing system Einstein@Home.Such searches are very computationally demanding, and would take decades to centuries on a single computer while still taking weeks or months on Einstein@Home.Thus, the search methods are specifically designed to ensure efficiency (Nieder et al. 2020).One key element is the use of constraints derived from optical observations.The companion's pulsar-facing side is heated by the pulsar wind, leading to a periodically varying optical light curve.This permits the orbital period P orb and other orbital parameters to be tightly constrained (for a feasible search the uncertainty ∆P orb needs to be less than a few milliseconds).In addition, because the sky position of the optical source is typically known to high precision (sub-milliarcsecond level), a search over position parameters is not needed.
Here we present the discovery and analysis of gamma-ray pulsations from PSR J1653−0158 in 4FGL J1653.6−0158.The pulsar is spinning very rapidly, at a rotational frequency of 508 Hz.The inferred surface magnetic-field strength is one of the lowest of all known pulsars.The discovery also confirms the 75 min orbital period.This very short orbital period raises interesting questions about the evolutionary path which created the system.
The paper is organized as follows.In Section 2, we describe the gamma-ray search, detection and analysis within LAT data.The optical analysis of the pulsar's companion, ra-dio pulsation searches, and a continuous gravitational-wave follow-up search are presented in Section 3. We discuss the results and conclude in Section 4.
Using the photon incidence angles and energies, we constructed a probability or weight for each photon, w j ∈ [0, 1], where j labels the photon: w j is the probability that the jth photon originated from the posited source, as opposed to a fore-or background source.These weights were computed by gtsrcprob, using the preliminary Fermi-LAT 8-year source catalog2 as a model for the flux within the RoI without performing a full spectral fit.Weighting the contribution of each photon to a detection statistic in this way greatly increases the search sensitivity (Kerr 2011), and the distribution of weights can be used to predict expected signal-tonoise ratios (Nieder et al. 2020).
The data set used here consisted of N = 354,009 photons, collected over a period of 3,542 days.The properties of the detection statistics (semicoherent power S 1 , coherent power P 1 , and H statistic) depend upon the lowest moments of the weights, which are These moments determine the ultimate sensitivity to a particular pulse profile and pulsed fraction, as given in Eq. ( 11) in Nieder et al. (2020).
Following the pulsar discovery, we extended this dataset to February 23, 2020 (MJD 58,902), using the latest P8R3_SOURCE_V2 IRFs (Bruel et al. 2018), a larger maximum zenith angle of 105 • , and using the Fermi-LAT Fourth Source Catalog (hereafter 4FGL, Abdollahi et al. 2020) as the RoI model for the photon probability weight computations.

Search
The binary-pulsar search methods are described by Nieder et al. (2020), which are a generalization and extension of the isolated-pulsar search methods from Pletsch & Clark (2014).
The sky position of the candidate optical counterpart is constrained to high precision in the Gaia catalog, so no astrometric search is required.The proper motion measured by Gaia for the optical counterpart was ignored for the search.

Orbital Constraints from Optical Observations
The orbital period estimate of Romani et al. (2014) was derived from Southern Astrophysical Research (SOAR), WIYN and Catalina Sky Survey (CSS) observations.These were augmented by new 350s SOAR Goodman High Throughput Spectrograph (GHTS) g , r , i exposures (63 g , 75 r , 42 i ) from MJD 56,516.184,and with the 300 s g , r and i exposures obtained by Kong et al. (2014) using the Wide Field camera (WFC) on the 2.5m Isaac Newton Telescope (INT) on La Palma.For these two data sets, the scatter about the light curve trends was appreciably larger than the very small statistical errors; we thus add 0.03 mag in quadrature to account for unmodeled fast variability and/or photometry systematics.To further refine the orbital period uncertainty, we obtained additional observations in u , g and i using the high-speed multi-band imager ULTRACAM (Dhillon et al. 2007) on the 4.2m William Herschel Telescope (WHT) on two nights (MJDs 57,170 and 57,195), covering six and three orbits of the binary system, respectively, with a series of 20 s exposures.Conditions were very poor on the first night with seeing > 5 arcsec, particularly at the beginning of the observation.We therefore only used the second night's data for the optical light curve modeling in Section 3.1, adding the latter half of the first night's observations for orbital period estimation.Finally, we obtained further INT+WFC exposures (23 g , 151 r , 45 i ) on MJD 57,988 -57,991.The g , r , i filter fluxes were referenced to in-field PanSTARRS catalog sources, and then converted to the SDSS scale.The u photometry was calibrated against an SDSS standard star observed on MJD 57,170.We estimate ∼ 0.05 mag systematic uncertainties in g , r and i , with uncertainties as large as ∼ 0.1 mag in u .
We constrained the orbital period using the multi-band Lomb Scargle periodogram method (VanderPlas & Ivezić 2015, excluding the u ULTRACAM data, as the modulation has very low signal-to-noise ratio in this band).To infer reasonable statistical uncertainties, we fit for and removed constant magnitude offsets, consistent with our estimated calibration uncertainties, between each night's observations in each band, and additionally rescaled the magnitude uncertainties to obtain a reduced chi-square of unity.This constrained the orbital period to P orb = 0.0519447518 ± 6.0 × 10 −9 days, where the quoted uncertainty is the 1σ statistical uncertainty.For the pulsation search, we chose to search the 3σ range around this value.
In Romani et al. (2014), the time of the pulsar's ascending node, T asc , was estimated from the photometric light curve.However, the optical maximum is distinctly asymmetric (see Section 3.1), which can bias orbital phase estimates.We therefore used the spectroscopic radial velocity measurements from Romani et al. (2014), folded at the orbital period obtained above, and fit the phase of a sinusoidal radial velocity curve, finding T asc = MJD 56513.47981± 2.1 × 10 −4 .However, as radial velocities may still be slightly biased by asymmetric heating, we elected to search a wide range around this value, corresponding to ±8σ.
For the projected semimajor-axis parameter x = a 1 sin i/c, we decided to start searching x ∈ [0, 0.1] s, with the intention to go to larger values in the case of no detection.For a pulsar mass of 1.6 M , this would cover the companion mass range up to 0.2 M and would include companion masses of all known "black-widow" systems as well as some of the lower mass "redback" systems (Roberts 2013;Strader et al. 2019).
Here, a 1 is the Die kpulsar's semimajor axis, i denotes the inclination angle, and c the speed of light.As described in Nieder et al. (2020), we expected x ∈ [0, 0.2] s based on the companion's velocity amplitude reported by Romani et al. (2014) and the masses expected for "spider" companions, i.e. black-widow or redback companions.

Search grids
To cover the relevant orbital-parameter space in {x, P orb , T asc }, we use optimized grids (Fehrmann & Pletsch 2014).These grids use as few points as possible still ensuring that a signal within the relevant space should be detected.Furthermore, they are able to cover the orbital-parameter space efficiently even though the required density depends on one of the orbital parameters, x.
Key to building an optimized grid is to know how the signal-to-noise ratio drops due to offsets from the true pulsar parameters.This is estimated using a distance metric on the orbital-parameter space (Nieder et al. 2020).In our case, the three-dimensional grid was designed to have a worst-case mismatch m = 0.2, i.e. not more than 20% of the (semicoherent or coherent) signal power should be lost due to orbitalparameter offsets.Of most relevance is that 99% of randomly injected orbital-parameter points have a mismatch below m = 0.04 to the closest grid point.
Due to the f -dependency of the required grid-point density, we search f in steps, and build the corresponding orbital grids prior to the start of the search on the computing cluster ATLAS in Hannover (Aulbert & Fehrmann 2008).

Searching
the 5-dimensional parameter space { f , ḟ , x, P orb , T asc } is a huge computational task with over 10 17 trials.Thus, the first (computing-intensive) search stages were performed on Einstein@Home, a distributed volunteer computing system (Allen et al. 2013).As done for radio pulsar searches previously, the search code utilizes the approximately 10,000 GPUs active on Einstein@Home for a computing speed-up of ∼ 10, comparing the runtimes on CPUs and GPUs.
The parameter space is divided into more than one million regions.Searching one of these is called a "work unit".These work units are sent to computers participating in Einstein@Home, and are searched when the computer is otherwise idle.Depending on the system, searching a work unit takes between half an hour and up to a few hours of computational time.In total, the search would have taken more than 50 years on a single computer, but using Einstein@Home it took less than two weeks.

Gamma-ray detection
The search process involves multiple stages in which semicoherent 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 ratio definitions, see Nieder et al. 2020).In the last semicoherent stage, a candidate found at a frequency of 1016 Hz had signal-to-noise ratio S 1 = 8.6, which we now associate with PSR J1653−0158.This was not the strongest candidate or far above the background of noise, but was among the ten most significant candidates in its work unit, and therefore passed on to the coherent stage.In the coherent stage, it was very significant, with signal-to-noise ratio P 1 /2 = 94.
The search follow-ups confirmed significant pulsations with period P ≈ 1.97 ms (or f ≈ 508 Hz), while the actual search revealed an alias at twice the pulsar frequency.This may be because the signal has significant power in the second harmonic.
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 by asymmetric heating (see Section 2.2.1).

Timing
The parameters used in the phase model to describe the pulsar's rotation are measured in a timing analysis.We use the timing methods as explained in Clark et al. (2017), which are an extension of the methods by Kerr et al. (2015).The basic principle is that the parameter space around the discovery parameters is explored using a Monte Carlo sampling algorithm with a template pulse profile.
To marginalize over the pulse profile template, we vary the template parameters as described in Nieder et al. (2019).In the case of PSR J1653−0158, we used a template consisting of two symmetrical, wrapped Gaussian peaks.We used constraints on the peaks' full-width at half maximum, such that the peaks must be broader than 5% of a rotation, and narrower than half a rotation.
Our timing solution over 11 years of LAT data is shown in Table 1.The folded gamma-ray data and the pulse profile are portrayed in Fig. 1.
The observed spin-down Ṗ is one of the lowest of all known pulsars.To estimate the intrinsic Ṗ we account for the Shklovskii effect (Shklovskii 1970), and the Galactic acceleration (see, e.g., Damour & Taylor 1991).The results are summarized in Table 1.The observed contribution due to the difference in Galactic acceleration of the Sun and the pulsar is computed with R Sun = 8.21 kpc, z Sun = 14 pc, and the Galactic potential model PJM17_best.Tpot (McMillan 2017), as implemented in their code 5 .For PSR J1653−0158, we used R J1653 = 7.48 kpc, and z J1653 = 367 pc, assuming d = 840 pc (see Table 2).The contributions parallel and perpendicular to the Galactic disk nearly cancel each other, so that the choice of the potential and its relevant parameters have a seemingly large effect on the actual small value of ṖGal , and can even change the sign.However, the overall kinematic contribution to the observed Ṗ is dominated by the Shklovskii term, and its uncertainty by the uncertainty in the distance estimate.The estimated intrinsic spin-down is Ṗint = 8.5 × 10 −22 s s −1 for distance d = 840 pc.

Optical Light Curve Modeling and System Masses
By modeling the optical light curves and radial velocities we can constrain the binary mass and distance and the system viewing angle.Comparing the individual filters between nights suggest small δm ≈ 0.05 shifts in zero points, consistent with the systematic estimates above.Correcting to match the individual filters, we then re-binned the light curve, placing the photometry on a regular grid with points spaced by δφ = 0.004, using the Python package Lightkurve; after excision of a few obviously discrepant points, we retain 248 u , 239 g , 220 r and 245 i 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 effect of hot spots on the companion surface, likely generated by precipitation of particles from the intrabinary shock (IBS) to companion magnetic poles (Sanchez & Romani 2017).All parameter values and errors are determined by MCMC modeling.The very shallow modulation of these light curves might normally be interpreted as indicating a small inclination i. However given the large companion radial velocity amplitude K = 666.9± 7.5 km s −1 , implying a mass function f (M) = 1.60 ± 0.05 M , measured by Romani et al. (2014), a small inclination would give an unphysical, large neutron star mass.As noted in that paper, the light curves and spectra show that a strong blue non-thermal veiling flux dominates at orbital minimum.With increasingly shallow modulation for the bluer colors, this is also evident in the present photometry.Thus the minimal model for this pulsar must include a non-thermal veiling flux.Although this is likely associated with the IBS, we model it here as a simple power law with form f ν = f A (ν/10 14 Hz) −p .This flux is nearly constant through the orbit, although there are hints of phase structure, e.g. in r and i at φ B = 0.72 (see Fig. 2).Any model without such a power law component is completely unacceptable.These fits prefer an A V slightly higher than, but consistent with the maximum in this direction (obtained by ∼ 300 pc; Green et al. 2019) 6 .
In Fig. 2, one notices that the orbital maximum is slightly delayed from φ B = 0.75, especially in the bluer colors.Such asymmetric heating is most easily modeled adding a polar hot spot with location (θ c , φ c ) and local temperature increase A c in a Gaussian pattern of width σ c ; when we include such a component, the fit improves greatly, with ∆χ 2 /DoF = −0.34.The Akaike Information Criterion comparison of the two models indicates that the model with a hot spot is preferred at the 10 −18 level, despite the extra degrees of freedom.We give the fit parameters for both models in Table 2.Note that with the fine structure near maximum, the model is not yet fully acceptable (χ 2 /DoF ∼ 1.4).More detailed models, including direct emission from the IBS or possibly the effects of companion global winds (Kandel & Romani 2020), may be needed to fully model the light curves.Such modeling would be greatly helped by light curves over an even broader spectral range, with IBS effects increasingly dominant in the UV, and low temperature companion emission better constrained in the IR.With many cycles we could also assess the reality (and stability) of the apparent fine structure and test for hot spot motion.
Our fit distance may be cross-checked with two other quantities.
(1) With the 4FGL energy flux f γ = 3.5 × 10 −11 erg cm −2 s −1 between 100 MeV and 100 GeV, our fit distance gives an isotropic gamma-ray luminosity L γ = 3 × 10 33 erg s −1 , in good agreement with the L γ ≈ (10 33 erg s −1 Ė) 1/2 heuristic luminosity law (Abdo et al. 2013), as a function of the spin-down power Ė.This luminosity is consistent with the model for direct radiative heating of the companion.(2) Our fit distance is also consistent with the model-independent, but lower accuracy, distance from the Gaia parallax.Thus, the 840 pc distance seems reliable, although systematic effects probably dominate over the rather small ∼ 50 pc statistical errors.
Armed with the fits, we can estimate the companion masses, correcting the observed radial velocity amplitude (fit with a K-star template) for the temperature-dependent weighting of the absorption lines across the companion face as in Kandel & Romani (2020).The results indicate substantial mass accretion, as expected for these ultra-short period systems.With the preferred Veiled+HS model the mass significantly exceeds 2.0 M , adding to the growing list of spider 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 on the heating pattern, can reduce the i uncertainty.

Radio pulsation searches
The pulsar position has been observed in radio multiple times.Several searches were performed before the gammaray pulsation discovery, and a few very sensitive follow-up searches afterwards.Despite the more than 20 observations with eight of the most sensitive radio telescopes, no radio pulsations have been found.
The results of the radio searches are given in Table 3. Observations are spread over 11 years, with observing frequencies ranging from 100 MHz up to 5 GHz.All orbital phases have been covered by most of the telescopes.Since there was no detection, the table also gives upper limits derived from the observations.For all but LOFAR, the data (both archival and recent) was folded with the gamma-ray-derived ephemeris, and searched only over dispersion measure.
The strictest upper limits on pulsed radio emission are 8 µJy at 1.4 GHz, and 20 µJy at 4.9 GHz.This is fainter than the threshold of 30 µJy that Abdo et al. (2013) use to define a pulsar to be "radio quiet".Note, that for the calculation of the limits we included the parts of the orbit where eclipses might be expected for spider pulsars.Thus, the limit constrains the maximum emission of the system, and not the maximum emission from the pulsar alone.

Continuous gravitational waves
We search for nearly monochromatic, continuous gravitational waves (GWs) from PSR J1653−0158, using data from the first7 and second8 observing runs of the Advanced LIGO detectors (The LIGO Scientific Collaboration et al. 2019).We assume that GWs are emitted at the first and second harmonic of the neutron star's rotational frequency, as would occur if the spin axis is misaligned with the principal axes of the moment of inertia tensor (Jones 2010(Jones , 2015)).
We employ two different analysis procedures, which yield consistent results.The first is frequentist, based on the multidetector maximum-likelihood F-statistic introduced by Cutler & Schutz (2005).The second is the Bayesian timedomain method (Dupuis & Woan 2005) as detailed by Pitkin et al. (2017), with triaxial non-aligned priors (Pitkin et al. 2015).Both methods coherently combine data from the two detectors, taking into account their antenna patterns and the GW polarization.The F-statistic search excludes data taken during times when the relevant frequency bands are excessively noisy.
The results are consistent with no GW emission.At twice the rotation frequency, the F-statistic 95% confidence upper limit on the intrinsic GW amplitude h 0 is 4.4 × 10 −26 .The 95%-credible interval upper limit from the Bayesian analysis on h 0 = 2C 22 is 3.0 × 10 −26 .At the rotation frequency (only checked with the Bayesian method) the 95% confidence upper limit on the amplitude C 21 is 6.6 × 10 −26 .
Since the dominant GW frequency might be mismatched from twice the rotation frequency (Abbott et al. 2019a), we performed an F-statistic search in a ±1 Hz band around this, with an extended ḟ -range.This yields larger upper limits on h 0 , with mean value of 1.3 × 10 −25 in 10 mHz-wide bands.Full details are given in the supplementary materials.
Our upper limits on h 0 at twice the rotation frequency may also be expressed as upper limits on the ellipticity of the pulsar (Abbott et al. 2019b).This is = 3.9 × 10 −8 × (h 0 /5 × 10 −26 ) × (10 45 g cm 3 /I zz ) × (840 pc/d), where I zz is the moment of inertia about the spin axis, and d is the distance.
As is case for most known pulsars, it is unlikely that our searches would have detected a GW signal.In fact, suppose that all of the rotational kinetic-energy losses associated with the intrinsic spin-down is via GW emission.Then assuming the canonical I zz = 10 45 g cm 3 , this would imply a "spin-down" ellipticity sd = 4.7 × 10 −10 , which is a factor ∼ 80 below our upper limit.

DISCUSSION & CONCLUSIONS
PSR J1653−0158 is the second binary pulsar (Pletsch et al. 2012), and the fourth MSP (Clark et al. 2018) to be discovered through periodicity searches of gamma rays.This pulsar is remarkable in many ways.It is only the second rotationally-powered MSP from which no radio pulsations have been detected.It is among the fastest-rotating known pulsars with spin frequency f = 508 Hz.The 75 min orbital period is shorter than for any other known rotation-powered pulsar, with the previous record being PSR J1311−3430 with a 93 min orbit (Pletsch et al. 2012).The inferred surface magnetic field is possibly the weakest, depending on the Shklovskii correction.
The discovery was enabled by constraints on the skyposition and orbital parameters from optical observations, together with efficient search techniques and the large computing power of the distributed volunteer computing system Einstein@Home.The detection proves that the optically variable candidate counterpart (Kong et al. 2014;Romani et al.Table 3. Summary of radio searches for PSR J1653−0158.The columns show the telescope used, the observed frequency range, the start time and data span, the range of orbital phases covered, the resulting limit on a pulsed component, and a reference with relevant details.NOTE-The orbital phase is given in orbits, and ranges > 1 indicate that more than one orbit has been observed.The considered maximum dispersion measure varies with the observing frequency from DM = 80 pc cm −3 at the lowest frequencies to DM = 350 pc cm −3 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.25 P, and a threshold signal-to-noise ratio S/Nmin = 7.

2014
) is indeed the black-widow-type binary companion to PSR J1653−0158, and it conclusively resolves the nature of the brightest remaining unidentified gamma-ray source, first found more than two decades ago (Hartman et al. 1999).
The distance to PSR J1653−0158, and its proper motion are well constrained.Gaia measurements of the parallax, ϖ = 1.88 ± 1.01 mas, imply a distance d = 530 +470 −200 pc.A consistent, but tighter constraint is given by our optical modeling with d = 840 +40 −40 pc.The proper motion (see Table 1) is also measured with good precision (Gaia and our timing are in agreement).PSR J1653−0158 has one of the lowest observed spinperiod derivatives of all known pulsars ( Ṗ = 2.4 × 10 −21 s s −1 ).The intrinsic Ṗ = 8.5×10 −22 s s −1 (accounting for Galactic ac-celeration and Shklovskii effects) is even smaller.In Fig. 3, PSR J1653−0158 is shown in a P-Ṗ diagram, alongside the known radio and gamma-ray pulsar population outside of globular clusters.
The intrinsic Ṗ can be used to estimate the pulsar's spin-down power Ė, surface magnetic-field strength B surf , magnetic-field strength at the light cylinder B LC , and characteristic age τ c .These are given in Table 1 for d = 840 pc.Constant lines of Ė, B surf , and τ c are displayed in Fig. 3 to show the distance-dependent ranges.
Spider pulsars in very-short-period orbits are difficult to discover with traditional radio searches.Even though we can now fold the radio data with the exact parameters, PSR J1653−0158 is still not visible.There are two simple explanations for the non-detection of radio pulsations.(1) PSR J1653−0158 is intrinsically radio-quiet, in that its radio beam does not cross the line-of-sight, or it has a very low luminosity.There is one other radio quiet MSP known (Clark et al. 2018).(2) Radio emission is blocked by material produced by the pulsar evaporating its companion.
The minimum average density of the companion 64 g cm −3 is very high, assuming a filled Roche lobe (Eggleton 1983).Using the filling factor from optical modeling, the average companion density 73 g cm −3 is even higher.The high density and the compact orbit suggest that the companion may be a helium white dwarf remnant, and that the system may have evolved from an ultracompact X-ray binary (Sengar et al. 2017;Kaplan et al. 2018).
The discovery of PSR J1653−0158 is the result of a multiwavelength campaign.The pulsar-like gamma-ray spectrum, and the non-detection of radio pulsations, motivated the search for a visible companion.This was subsequently discovered in optical and X-ray observations.Further optical observations provided constraints on the orbital parameters which were precise enough to enable a successful gammaray pulsation search.

ACKNOWLEDGMENTS
We are deeply grateful to the thousands of volunteers who donated their computing time to Einstein@Home, and to those whose computers first detected PSR J1653−0158: Yi-Sheng Wu of Taoyuan, Taiwan; and Daniel Scott of Ankeny, Iowa, USA.
This work was supported by the Max-Planck-Gesellschaft (MPG), by the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Grant, No. PL 710/1-1 (PI: Holger J.The ULTRACAM photometry was obtained as part of program WHT/2015A/35.The William Herschel Telescope is operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.Based on observations made with the Isaac Newton Telescope (program I17BN005) operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.This paper makes use of data obtained from the Isaac Newton Group of Telescopes Archive which is maintained as part of the CASU Astronomical Data Centre at the Institute of Astronomy, Cambridge.
We acknowledge support of the Department of Atomic Energy, Government of India, under project no.12-R&D-TFR-5.02-0700for the GMRT observations.The GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research, India.The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS).We acknowledge financial support from the "Programme National Hautes Energies" (PNHE) of CNRS/INSU, France.This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT) under project code LC7_018.LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. FAST is a Chinese national mega-science facility, built and operated by NAOC.Partly based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg.
The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis.These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.
Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.This work performed in part under DOE Contract DE-AC02-76SF00515.
The authors thank the LIGO Scientific Collaboration for access to the data and gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, and the Max-Planck-Society (MPS) for support of the construction of Advanced LIGO.Additional support for Advanced LIGO was provided by the Australian Research Council.This research has made use of data, software and/or web tools obtained from the LIGO Open Science Center (https://losc.ligo.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration, to which the authors have also contributed.LIGO is funded by the U.S. National Science Foundation.Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.

Figure 1 .
Figure 1.Integrated pulse profile and phase-time diagram of PSR J1653−0158, showing two identical rotations.Top: The histogram shows the weighted counts for 50 bins.The orange curve indicates the pulse-profile template with the highest signal power, and the transparent black curves represent 100 templates randomly selected from the Monte Carlo samples after the chain stabilized, to indicate the uncertainty on the profile.The dashed blue line denotes the source background.Bottom: Each point represents the pulsar's rotational phase at emission of a photon, with the intensity indicating the photon's probability weight.Note that PSR J1653−0158 received more exposure between MJDs 56,600 and 57,000 when the LAT pointed more often towards the Galactic center.

Figure 2
Figure 2. u , g , r , and i light curves for PSR J1653−0158, with the best-fit model curves.Note the flat minima and decreasing modulation for bluer colors, a consequence of the hard spectrum veiling flux.Two identical cycles are shown for clarity.

Figure 3 .
Figure3.Newly detected PSR J1653−0158 on a P-Ṗ diagram of the pulsar population outside of globular clusters.The MSP population is shown magnified in the inset.LAT pulsars are marked in green (isolated by a cross and binary by a circle).Non-LAT pulsars in the ATNF are marked in gray (isolated by a plus and binary by a square).The lines show constant surface magneticfield strength (dashed-dotted), characteristic age (dotted), and spindown power (dashed).The spin period and intrinsic spin-period derivative of PSR J1653−0158 are marked by the orange star.The transparent stars indicate the (distance-dependent) maximum and minimum intrinsic spin-period derivative according to the distance estimated from our optical models.

Table 2 .
Light curve fit results for PSR J1653−0158.Parameters from the best-fit light curve/radial velocity models, with and without a surface hot spot, including MCMC errors.
Pletsch) and by National Science Foundation grants 1104902 and 1816904.L.N. was supported by an STSM Grant from COST Action CA16214.C.J.C. and R.P.B. acknowledge support from the ERC under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 715051; Spiders).V.S.D. and ULTRACAM are supported by the STFC.R.W.R. and D.K. were supported in part by NASA grant 80NSSC17K0024.S.M.R. is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center award 1430284 and the NASA Fermi GO Award NNX16AR55G.Fermi research at NRL is funded by NASA.J.W.T.H. is an NWO Vici fellow.