PSR J0952-0607: The Fastest and Heaviest Known Galactic Neutron Star

We describe Keck-telescope spectrophotometry and imaging of the companion of the ``black widow"pulsar PSR~J0952$-$0607, the fastest known spinning neutron star (NS) in the disk of the Milky Way. The companion is very faint at minimum brightness, presenting observational challenges, but we have measured multicolor light curves and obtained radial velocities over the illuminated ``day"half of the orbit. The model fits indicate system inclination $i=59.8\pm 1.9^\circ$ and a pulsar mass $M_{NS} = 2.35\pm 0.17 M_\odot$, the largest well-measured mass found to date. Modeling uncertainties are small, since the heating is not extreme; the companion lies well within its Roche lobe and a simple direct-heating model provides the best fit. If the NS started at a typical pulsar birth mass, nearly $1 M_\odot$ has been accreted; this may be connected with the especially low intrinsic dipole surface field, estimated at $6\times 10^7$G. Joined with reanalysis of other black widow and redback pulsars, we find that the minimum value for the maximum NS mass is $M_{\rm max}>2.19 M_\odot$$(2.09 M_\odot)$ at $1\sigma$$(3\sigma)$ confidence. This is $\sim 0.15 M_\odot$ heavier than the lower limit on $M_{\rm max}$ implied by the white-dwarf--pulsar binaries measured via radio Shapiro-delay techniques.


INTRODUCTION
Pulsar PSR J0952−0607 (hereafter J0952) was discovered by Bassa et al. (2017) with a spin period of P s = 1.41 ms, making it the fastest-spinning pulsar in the disk of the Milky Way. It is a "black widow" (BW) pulsar with a low-mass (substellar) companion being irradiated and evaporated by the pulsar luminosity. Nieder et al. (2019) subsequently detected it as a gamma-ray pulsar, and obtained additional radio timing and optical photometry that allowed initial fits for the system properties. In particular, they measure a low period derivativė P obs = 4.6 × 10 −21 s s −1 , which places an upper limit on the surface dipole field of 8.2 × 10 7 G, among the 15 lowest-known pulsar magnetic fields (B), even before the Shklovskii (1970) correction. Since in addition their optical photometry suggests a large (> 5 kpc) distance, and timing data gave a best-fit (albeit low-significance) P =Ṗ obs − 2.43 × 10 −21 P µ 2 mas/yr d kpc , may be substantial, reducingṖ to 2.0 × 10 −21 s s −1 and thus B to B i ≈ 6.1 × 10 7 G. Not all pulsars have a proper motion allowing an intrinsic B i estimate, but only ten pulsars have a lower B i in the ATNF catalog (Manchester et al. 2005, http://www.atnf.csiro.au/research/pulsar/psrcat). The measurement of high neutron star (NS) masses in some white dwarf (WD)-NS millisecond pulsar binaries, via radio-measured Shapiro delay, has been very influential in driving thinking about the dense-matter equation of state. Two systems have best-estimate masses barely exceeding 2.0 M : J0348+0432 at 2.01 ± 0.04 M (Antoniadis et al. 2013) and PSR J0740+6620 at 2.08 ± 0.07 M (Fonseca et al. 2021), both 1σ uncertainties. However, the maximum mass that an NS can reach must depend on its binary evolutionary history. For example, inspection of the double NS-NS binaries alone would lead one to conclude that the typical NS mass is arXiv:2207.05124v1 [astro-ph.HE] 11 Jul 2022 ∼ 1.35 M , with a maximum of ∼ 1.45 M ; the millisecond pulsars in these systems are only mildly recycled, with spin periods of tens of ms and moderate magnetic fields. One must therefore examine a variety of pulsar binary classes in a quest to find the most massive NSs.
It has long been argued (e.g., Bisnovatyi-Kogan & Komberg 1974;Romani 1990) that binary accretion may reduce young pulsar terraGauss fields to millisecondpulsar values. While the mechanism is unclear, ranging from simple burial to heating-driven ohmic decay (see Mukherjee 2017, for a review), it seems likely that the amount of mass accretion or its duration play a role in controlling the degree of field reduction. Thus, it is interesting to measure the masses of millisecond pulsars, and the fastest-spinning lowest-field pulsars are good candidates for substantial accretion and high NS mass. Since evolution is driven by angular momentum loss and irradiation bloating, rather than by nuclear evolution, long periods of sub-Eddington rate accretion are expected for the progenitors of "spider" [BW and "redback" (RB) millisecond pulsar] binaries before the start of the pulsar-driven evaporation phase; these shortperiod binaries may host very massive NSs. Thus, J0952 is a particularly attractive candidates for further investigation.
However, Nieder et al. (2019) state, "Unfortunately, the optical counterpart of PSR J0952−0607 is too faint (r ≈ 23 at quadrature when the radial velocity is highest) for spectroscopic radial velocity measurements to be feasible even with 10 m class telescopes." This is nearly true, but the exceptional P s and B inspired our intensive campaign of J0952 Keck LRIS imaging and spectroscopy. Modeling these data, we find that the NS mass is indeed the largest well-measured value to date. Combining this mass with those from new modeling of other spider binaries, we show that these objects put a strong lower bound on the maximum NS mass, appreciably above the values inferred from radio Shapiro delay measurements. This should have deep implications for the dense-matter equation of state.

OBSERVATIONS
With its P b = 6.42 hr orbital period and faint magnitude, the effective study of J0952 requires substantial amounts of large-telescope time, ideally with superior seeing. Our campaign uses the Keck-I 10 m telescope and Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995). We report here on two dual-color (typically g/i with a few r/i pairs) imaging runs and six spectroscopy campaigns (Table 1). For the latter we used the 5600Å dichroic, the 600 l/4000Å blue grism, and the 400 l/8500Å red grating, covering ∼ 3300-10,500Å with dispersions of 0.63Å pixel −1 (blue side) and 1.2Å pixel −1 (red side). Typical exposures were 600-900 s. All spectra were processed with the LPipe reduction software (Perley 2019). The LRIS red arm thick chip collects many cosmic rays, so additional editing was needed to clean the red-arm spectra. All companion and comparison star spectra are available as a DbF (Data behind Figures) tarball.
With the atmospheric dispersion corrector (ADC) being used, we could rotate the 1 -wide slit away from the parallactic angle (Filippenko 1982) to simultaneously measure a nearby brighter G1 star with known Pan-STARRS2 (PS2) magnitudes. This allows us to monitor the system throughput and wavelength solution among frames. In particular, since this PS2 star has known and stable magnitudes, we integrate the spectra over the SDSS standard ugriz filter bands using the sbands IRAF script, and calibrate with the comparison star's catalog magnitudes (converted to the SDSS system using the prescription of Finkbeiner et al. 2016) to obtain light curves with absolute fluxes, up to small gray shifts from differential slit losses.
Our spectroscopic strategy was to observe J0952 for half-nights ideally covering optical maximum brightness. To cover the critical phases at low airmass under dark moon conditions required a careful selection of observing nights. Unfortunately, none of the nights had particularly good seeing; indeed, several additional halfnights requested for this program were completely lost to weather and COVID-19 shutdowns. Nevertheless, the night-to-night consistency in the flux scale and radialvelocity (RV) scale, aided by the comparison star, has let us make combined fits of all the data, allowing good model constraints.
As noted, imperfect position angle (PA) and slit placement can give small differential slit losses. These we correct via gray shifts to the sbands magnitudes, matching the imaging fluxes in the overlap regions and between epochs. These ≤ −0.2 mag shifts increase the companion flux; this is sensible, as one tends to ensure that the (brighter) comparison star is well covered by the slit.  In the phase window φ B = 0.1-0.4 the companion is too faint to determine spectrophotometric colors. We therefore collected LRIS dual-band (g/i with a few g/r) direct images covering minimum brightness. We established a reference grid of PS2 catalog stars, with po- Figure 3. RV measurements for J0952. Dot size is proportional to the correlation coefficient R; low-significance correlations with R < 2 are marked with black crosses and are excluded from all fits. Three higher-significance points marked with green crosses appear as outliers and are optionally excluded for the "trimmed" fit. The background curve is the best-fit RV curve for the heating model fit to the light curve in Figure 1. The lower panel shows the fit residuals.
sitions set relative to the pulsar companion, measured near maximum brightness, and extracted forced aperture photometry at these fixed positions. PS2 magnitudes were converted to the SDSS scale and used to determine the companion magnitude. Direct photometric errors were combined in quadrature with zero-point systematic errors determined from scatter in the catalog star residuals; near minimum the photometric errors strongly dominate the final uncertainty. The resulting light curve covers a large (> 5 mag) range, an observational challenge. Image artifacts seem to produce a few modest outliers, and companion flare activity may occasionally brighten the source. We discuss trimming such points below.
We measure the RV amplitude by IRAF rvsao (Kurtz & Mink 1998) cross-correlation with a G1 template, using the 3800-9500Å range, excluding wavelengths near strong telluric features. The spectra were resampled to 8192 wavelengths and filtered in Fourier space with a lower cosine-bell filter cutoff from 20-100 and the upper cutoff running from 1000-1800. The average companion spectrum near φ B = 0.75 is shown in Figure 2. The pulsar and comparison star were correlated using identical parameters. For the first two epochs (2018B and 2019A), a different slit PA led to the use of a cooler (∼ K5) comparison star. A small velocity offset was found between these data and the G1 comparison-star correlations used for the other epochs. In practice, the red K5 star produced poorer RV solutions, and lacked sufficient flux for good g and u calibration; thus, the 2018B and 2019A data do not contribute spectrophotometric points to Figure 1. The G1 comparison-star velocity was always measured to σ v < 5 km s −1 , but pulsar companion velocity uncertainties varied greatly from as little as 10 km s −1 near "midday" (φ B = 0.75) to 30-50 km s −1 toward the "dawn' and "dusk' phases. In addition, spurious correlation peaks appeared well outside the plausible companion-velocity range for many lower signal-to-noise-ratio (S/N) spectra.
Thus, we measured in two passes, first correlating all spectra over a broad velocity range, and fitting the results with a simple sinusoid RV of amplitude K and offset Γ at the radio ephemeris orbital phasing. We then correlated again, accepting only results within ±100 km s −1 of the preliminary fit. This greatly decreased the number of outlier points, and allowed us to recover several weaker RV measurements near the plausible companion solution. For each spectrum the velocity was referenced to that of the comparison star, which varied between exposures, evidently driven by details of the pipeline wavelength solution and extraction aperture. We confirm that these variations are in the pulsar, as well; the scatter in the final RVs increases significantly if the points are not referenced to the comparison star. Small overall RV zero-point shifts were also found in several epochs, again likely owing to differences in the apertures. Figure 3 shows all RV values with a finite correlation coefficient (R > 0). The point size is scaled to R; a number of obvious outliers with very small R remain.

PHOTOMETRIC/RADIAL-VELOCITY FITTING
Our light-curve fits are performed with an outgrowth of the ICARUS light-curve model (Breton et al. 2012) with the additions described by Kandel & Romani (2020) and Romani et al. (2021); the main parameters of interest are the orbital inclination i, the companion Roche lobe fill factor f 1 , the heating luminosity L H , the companion's night-side temperature T N , and the system distance d. In addition we can fit for the extinction A V to the source and a set of "bandpass calibrations," for possible imperfections in the photometric zero points.
The MCMC fits and Bayseian parameter estimation are performed by Multinest sampling (Feroz et al. 2009) using its Python implementation pymultinest (Buchner 2016). In addition to the photometry, the fits use the pulsar kinematic parameters of Nieder et al. (2019). The light-curve fit is also (weakly) dependant on the companion center of mass RV K CoM . We do not fit for this with ICARUS; instead, we initially fit light curves using the observed spectroscopic K and then refit with K CoM determined from the spectroscopic fit (below), iterating to convergence.
The fit χ 2 is increased by a handful of measurements several σ from the models. The synthesized u photometry has several such points owing to the very low S/N (per resolution element) of the spectra in this band. The most important outliers are a few bright g and i direct-imaging points during the "night" phase. Many spider binaries exhibit occasional short-term (∼ 1000 s) light-curve flares. Since the quiescence light curve is needed for the light-curve fits, we optionally excise the few points well above the model curves (marked in Figure 1). Table 2 gives the fit values with these points first excluded, and then included. Note that inclusion slightly decreases the fit i and increases the mass (by ∼ 0.7σ), while giving much larger χ 2 . We therefore adopt the fit to the "trimmed" dataset as conservative.
Extinction in the direction of J0952 is estimated from the three-dimensional dust map of Green et al. (2018), reaching its maximal E(g − r) = 0.05 +0.03 −0.02 mag (A V = 0.14 +0.08 −0.05 mag) by 0.3 kpc. Treating A V as a free parameter in the MCMC fits, we find 0.16 ± 0.04 mag, consistent with the dust-map value. Since we rely on the PS2 griz catalog magnitudes of our comparison star, there are ∼0.02 mag zero point uncertainties (lacking PS2 u, we rely on the extrapolation of Finkbeiner et al. 2016). Leaving the passband calibrations free in the ICARUS fit to the trimmed dataset, we find ∆u = +0.07 +0.08 −0.10 , ∆g = +0.05 +0.05 −0.06 , ∆r = +0.02 ± 0.04, ∆i = −0.02 ± 0.04 and ∆z = −0.02 ± 0.04 mag, all consistent with zero. We therefore believe that our initial calibration was quite good. For the untrimmed data, the outliers produce a large χ 2 /DoF and drive these color terms to larger values, notably ∆ugriz = +0.27, +0.23, +0.03, −0.12, −0.15 mag. While including these passband shifts does not strongly change either χ 2 or the fit parameters, it does slightly increase the MCMC uncertainty ranges of other parameters (most notably d and L H ); we thus retain these in the fit, although we omit them from Figure 4.
The fits all find i ≤ 60 • , with f 1 = 0.79 (companion underfilling its Roche lobe), T N ≈ 3000 K, day side heated to ∼ 6200 K, and d ≈ 6 kpc, larger than the dispersion-measure estimates. These values are very similar to those of Nieder et al. (2019), but with substantially higher precision. Interestingly, simple directheating models always provide the best statistical fit; adding additional effects required for many other spider binaries, such as companion hot spots and surface winds, does not significantly improve the fit. This means that our results are particularly immune to model-dependant systematics, and the fit-parameter errors (Table 2) are dominated by measurement uncertainties. We next compare the Keck spectroscopy with model RVs computed from photometric-fit companion models. As emphasized by Linares et al. (2018) and discussed by KR20, different species' temperature sensitivities cause varying line equivalent widths (EWs) across the face of the companion. As the model parameters change, the  surface heating changes and the line EWs vary. For J0952, metal lines dominate the cross-correlation signal (as appropriate for our G1 correlation spectrum); we apply metal line equivalent width EW(T ) weighting to correct center-of-mass RVs to observed center-of-light values (see KR20). If we (inappropriately) apply Balmer EW(T ) weighting, the RV amplitude increases to 384 ± 5km s −1 , so again our choice is conservative.
We always exclude the low-significance correlations with R < 2 in the RV MCMC fits (marked with black crosses in Figure 3; note that some crossed points follow the expected curve, but many are quite discrepant). The remaining χ 2 is in fact dominated by three outlier points with small uncertainties (marked with green crosses in Fig. 3). In our "trimmed" fit for the mass we exclude these three points. Again, this is conservative, as with them the K increases slightly (by 0.4σ). The lower section of Table 2 gives the RV fit parameters and their uncertainties. The mass difference in the "All" column is dominated by the smaller i of the "All" photometric fit.

MASS IMPLICATIONS AND CONCLUSIONS
In the quest for precision mass measurements from spider binaries, J0952 has the advantage that the complex modeling effects required to fit other spiders (hot spots, wind flows, extreme gravity darkening) are not needed. This is a natural consequence of its relatively large P b , weak pulsar heating, and low Roche-lobe fill factor. On the other hand, the large distance, weak heating, and low f 1 mean that the companion is quite faint, barely bright enough for measurement with present 10 m telescope systems. Our intensive Keck campaign has managed to measure the binary properties, constraining the masses with sufficient precision that the lower limit on M NS contributes to equation-of-state constraints.
We can examine the minimum required M max by plotting the Gaussian probability distribution functions with M i , σ i for observed heavy NSs in Figure 5. Here we compare the latest measurements of the > 2 M "radio selected" WD-NS MSP J0348+0432 and J0740+6620 with the measurements extracted from light-curve and RV measurements of spider pulsar (BW and RB) companions. For the latter we use the most conservative (best-fit) models with all heating effects as summarized by Kandel & Romani (2022). To estimate the minimum value of M max required by these observations, we assume that the NS mass distribution is flat M 0 above M 1 = 1.8 M to some cutoff value (for this > 2 M sample, results are insensitive to the underlying distribution; values change by < 0.01 M for M −1 ). The log(likelihood) for n measurements is These curves are plotted in the lower panel of Figure  5. We list the medians and lower bounds at various confidence levels on M max , from the likelihood ratio, for the different sample sets.
Since the spiders are appreciably heavier than the WD-NS binaries, the 1σ M max bound increases by ∼ 0.13 M , or ∼ 0.17 M with J0952. We have checked these bounds with several different estimators. For example a bias-corrected, accelerated bootstrap analysis (kindly supplied by B. Efron) of the seven masses gives M max = 2.34 M with a 2σ lower bound of 2.17 M . Numerical simulations with random draws and our observed errors give similar values. Of course, J0952 would constrain M max even more strongly if its mass uncertainty could be reduced; with σ = 0.05M , we would infer M max > 2.30 M (1σ) and > 2.20 M (3σ).
From studies of the millisecond pulsar population, it seems that the initial masses M i of these pulsars are bimodal, with a dominant component of M i ≤ 1.4 M and a 20% subpopulation with M i ≈ 1.8 M (Antoniadis et al. 2016, and references therein). This complicates use of the present mass to constrain mass accretion and test evolutionary models (e.g., of magnetic field reduction). For J0952 at our 1σ lower limit we infer that at least 0.5 M and more likely ∼ 1 M has been accreted. Assuming a start from 1.4 M , three of our BWs have particularly high accreted mass (J0952, ∆M ≈ 0.95 M ; J1653−0158, ∆M ≈ 0.77 M ; J1810+1740, ∆M ≈ 0.73 M ). After the Shklovskii (1970) corrections, these have some of the lowest pulsar dipole fields known (estimated as 6.1×10 7 G, 3.9×10 7 G, and 7.7×10 7 G, respectively). More M and B i measurements, as well as additional analysis, are needed to see if this is coincidental or causal.
Finally, we note that with a central value of 2.35 M , J0952 provides the most severe constraints on the densematter equation of state. This remains true even considering the lower bound, and this bound can be adopted with high confidence given the strong preference for a simple heating model. Of course, we would like an even tighter mass measurement of this especially important system. Improved photometry in blue filters at optical minimum may be feasible with 10 m-class telescopes and excellent conditions, but improved RVs likely await the 30 m telescope era.
We thank B. Efron for a discussion of bounds from a data sample and for the bootstrap analysis, and the anonymous referee whose probing questions led to some improvements in the paper. The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and NASA; it was made possible by the generous financial support of the W. M. Keck Foundation. We are grateful for the excellent assistance of the observa-tory staff, as well as the Keck Time Allocation Committee and Keck scheduler for accommodating our requests for specific half-nights. D.K. and R.W.R. were supported in part by NASA grants 80NSSC17K0024 and 80NSSC17K0502. A.V.F.'s group received generous financial assistance from the Christopher R. Redlich Fund, the TABASGO Foundation, and the U.C. Berkeley Miller Institute for Basic Research in Science (where A.V.F. was a Miller Senior Fellow).