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

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 derivative ${\dot{P}}_{\mathrm{obs}}=4.6\times {10}^{-21}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, which places an upper limit on the surface dipole field of 8.2 × 10⁷ 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) proper motion of ∼10 mas yr⁻¹, the correction to find the intrinsic $\dot{P}$,

$$\begin{eqnarray}&&\dot{P}={\dot{P}}_{\mathrm{obs}}-2.43\times {10}^{-21}P\,{\mu }_{\mathrm{mas}/\mathrm{yr}}^{2}{d}_{\mathrm{kpc}},\end{eqnarray} \tag{ 1 }$$

may be substantial, reducing $\dot{P}$ to 2.0 × 10⁻²¹ s s⁻¹ and thus B to B_i ≈ 6.1 × 10⁷ G. Not all pulsars have a proper motion allowing an intrinsic B_i estimate, but only 10 pulsars have a lower B_i in the ATNF catalog ³ (Manchester et al. 2005).

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 ∼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 milliseconds 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 millisecond-pulsar 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 short-period binaries may host very massive NSs.

Thus, J0952 is a particularly attractive candidate 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.

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 the 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–10500 Å with dispersions of 0.63 Å pixel⁻¹ (blue side) and 1.2 Å pixel⁻¹ (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 in a TAR archive as the data behind Figure 1.

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 Sloan Digital Sky Survey (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 half-nights requested for this program were completely lost to weather and COVID-19 shutdowns. Nevertheless, the night-to-night consistency in the flux scale and radial-velocity (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 positions 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 (Figure 1) 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.

Zoom In Zoom Out Reset image size

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 to 100 and the upper cutoff running from 1000 to 1800. The average companion spectrum near ϕ_B = 0.75 is shown in Figure 2. The pulsar and companion 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⁻¹, but pulsar companion velocity uncertainties varied greatly from as little as 10 km s⁻¹ near "midday" (ϕ_B = 0.75) to 30–50 km s⁻¹ 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.

Zoom In Zoom Out Reset image size

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⁻¹ 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.

Zoom In Zoom Out Reset image size

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 filling factor f₁, the heating luminosity L_H , the companion's nightside 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 Markov Chain Monte Carlo (MCMC) fits and Bayesian 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) dependent 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 χ² 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 χ². We therefore adopt the fit to the "trimmed" data set as conservative.

Table 2. Light-curve/RV-fit Results for J0952 ^a

Parameters	Trimmed	All
i ( deg)	${59.8}_{-1.9}^{+2.0}$	${58.5}_{-1.8}^{+1.9}$
f1	0.79 ± 0.01	0.77 ± 0.01
LH /1034 (erg s−1)	${3.81}_{-0.43}^{+0.46}$	${6.22}_{-0.77}^{+0.88}$
TN(K)	${3085}_{-80}^{+85}$	${3206}_{-95}^{+100}$
dkpc	${6.26}_{-0.40}^{+0.36}$	${7.60}_{-0.82}^{+0.74}$
χ2/DoF	286/(298-11)[1.00]	451/(314-11)[1.49]
KCoM (km s−1)	376.1 ± 5.1	379.1 ± 6.8
MNS(M⊙)	2.35 ± 0.17	2.50 ± 0.20
MC(M⊙)	0.032 ± 0.002	0.034 ± 0.002
χ2/DoF	55/(40-2)[1.4]	90/(43-2)[2.2]

Note.

^a Also fitted: A_V , Δu, Δg, Δr, Δi, and Δz. See the text.

Download table as:  ASCIITypeset image

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.02}^{+0.03}$ mag (${A}_{V}={0.14}_{-0.05}^{+0.08}$ 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 data set, we find ${\rm{\Delta }}u=+{0.07}_{-0.10}^{+0.08}$, ${\rm{\Delta }}g=+{0.05}_{-0.06}^{+0.05}$, Δ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 χ²/degrees of freedom (DoF) and drive these color terms to larger values, notably Δugriz = +0.27, +0.23, +0.03, −0.12, and −0.15 mag. While including these passband shifts does not strongly change either χ² 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.

Zoom In Zoom Out Reset image size

The fits all find i ≤ 60°, with f₁ = 0.79 (companion underfilling its Roche lobe), T_N ≈ 3000 K, dayside 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 direct-heating 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 ± 5 km s⁻¹, 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 χ² is in fact dominated by three outlier points with small uncertainties (marked with green crosses in Figure 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.

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 filling factor. On the other hand, the large distance, weak heating, and low f₁ 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 D. Kandel & R. W. Romani (2022, in preparation). To estimate the minimum value of ${M}_{\max }$ required by these observations, we assume that the NS mass distribution is flat M⁰ above M₁ = 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⁻¹). The log(likelihood) for n measurements is

$$\begin{eqnarray*}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & = & -n\,\mathrm{ln}({M}_{\max }-{M}_{1})\\ & & +{{\rm{\Sigma }}}_{i=1}^{n}\mathrm{ln}\left[\mathrm{erf}\left(\displaystyle \frac{{M}_{i}-{M}_{1}}{{2}^{1/2}{\sigma }_{i}}\right)-\mathrm{erf}\left(\displaystyle \frac{{M}_{i}-{M}_{\max }}{{2}^{1/2}{\sigma }_{i}}\right)\right].\end{array}\end{eqnarray*}$$

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.

Zoom In Zoom Out Reset image size

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}_{\odot }$ 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.05 M_⊙, we would infer ${M}_{\max }\gt 2.30\,{M}_{\odot }(1\sigma )$ 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⁷ G, 3.9 × 10⁷ G, and 7.7 × 10⁷ 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 dense-matter 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 observatory 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 80NSSC21K0896. 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).