XMM-Newton Detection and Spectrum of the Second Fastest Spinning Pulsar PSR J0952−0607

Despite sensitivity to and searches for pulsars with even higher spin rates, the fastest pulsars known to date among nonaccreting pulsars are PSR J1748−2446ad (with a spin rate ${\nu }_{{\rm{s}}}=716$ Hz; Hessels et al. 2006), PSR J0952−0607 (${\nu }_{{\rm{s}}}=707\,\mathrm{Hz};$ Bassa et al. 2017), and PSR B1937+21 (${\nu }_{{\rm{s}}}=642\,\mathrm{Hz};$ Backer et al. 1982) and among accreting pulsars are 4U 1608−52 (${\nu }_{{\rm{s}}}\,=620$ Hz), SAX J1750.8−2900 (${\nu }_{{\rm{s}}}=601$ Hz), and IGR 00291+5934 (${\nu }_{{\rm{s}}}=599$ Hz). The observed spin rates are well below the theoretical limit of ∼2000 Hz (Cook et al. 1994; Haensel et al. 1999). This suggests a mechanism that prevents fast rotation (Chakrabarty et al. 2003; Chakrabarty 2008; Papitto et al. 2014; Patruno et al. 2017; Gittins & Andersson 2019), such as mechanisms associated with gravitational wave emission since their torques depend strongly on spin rate.

One particular mechanism that is of great interest is that associated with the r-mode fluid oscillation, because it can be a strong source of gravitational waves (Andersson 1998; Friedman & Morsink 1998) via the Chandrasekhar–Friedman–Schutz instability (Chandrasekhar 1970; Friedman & Schutz 1978; see also Chugunov 2017). The strength of the r-mode instability is characterized by the balance between the timescale of mode growth by gravitational wave emission ${t}_{\mathrm{gw}}$ (which depends strongly on spin frequency, i.e., ${t}_{\mathrm{gw}}\propto {\nu }_{{\rm{s}}}^{-6}$) and timescale of viscous damping (which is temperature-dependent; e.g., ${t}_{\mathrm{visc}}\propto {T}^{2};$ Andersson & Kokkotas 2001). However, our theoretical picture of r-modes is severely problematic (see, e.g., Ho et al. 2011; Haskell et al. 2012; Chugunov et al. 2017). Thus it is vital to identify fast and hot neutron stars, which can be used to better understand the nature of the r-mode mechanism. Here we report the X-ray detection of the second fastest pulsar, PSR J0952−0607, which provides possibly the strongest r-mode constraints for millisecond pulsars.

PSR J0952−0607 was discovered at radio frequencies using LOFAR during a targeted search of gamma-ray sources detected by Fermi but not associated with other known sources (Bassa et al. 2017). Optical observations identify the binary companion of the pulsar, and the companion star's low mass and short 6.42 hr orbital period suggest PSR J0952−0607 is in a black widow system, where the pulsar wind irradiates and evaporates the companion star. The position of PSR J0952−0607 used for radio timing is that determined from its optical counterpart i.e., (R.A., decl.[J2000]) =(09^h52^m08319, $-06^\circ 07^{\prime} 23\buildrel{\prime\prime}\over{.} 49$). The distance is determined from the measured dispersion measure ($\mathrm{DM}=22.4\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$) and found to be either d = 0.97 kpc or 1.74 kpc, depending on which model of Galactic electron distribution is used (NE2001 or YMW16, respectively). Henceforth, we assume a distance of 1.74 kpc, unless otherwise noted. Bassa et al. (2017) estimate an interstellar absorption column ${N}_{{\rm{H}}}=4\times {10}^{20}$ cm⁻² from a Galactic extinction model and distance, which agrees with ${N}_{{\rm{H}}}=3.9\times {10}^{20}$ cm⁻² estimated from H i in the direction of PSR J0952−0607 (Dickey & Lockman 1990). We estimate a somewhat larger ${N}_{{\rm{H}}}\,={6.7}_{-2.0}^{+2.9}\times {10}^{20}$ cm⁻² (90% confidence level) using the empirical relation between ${N}_{{\rm{H}}}$ and DM from He et al. (2013). Since PSR J0952−0607 is in a short orbital period black widow system, a higher ${N}_{{\rm{H}}}$ could be measured, e.g., as might be possible using dispersion measure variations during radio eclipses such as that reported by Main et al. (2018) for PSR B1957+20. Using a short 4.6 ks exposure with Swift-XRT, Bassa et al. (2017) obtain a 3σ upper limit on the 0.3–10 keV flux f_0.3−10 < 1.1 × 10⁻¹³ erg cm⁻² s⁻¹, which corresponds to a X-ray luminosity limit of L < 1.1 × 10³¹ erg s⁻¹ at 0.97 kpc or 3.6 × 10³¹ erg s⁻¹ at 1.74 kpc.

In Section 2, we describe the XMM-Newton observation of PSR J0952−0607 and our procedure for processing the data. Sections 3 and 4 give details of our spectral fitting and variability search analyses. In Section 5, we summarize and discuss our results and use them to place constraints on the r-mode instability.

XMM-Newton observed PSR J0952−0607 on 2018 May 4 (ObsID 0821520101) for 71.2 ks with EPIC in full frame imaging mode and using the thin optical filter. Figure 1 shows the MOS2 and pn images of the field around PSR J0952−0607. We process this data using SAS 17.0.0 (Gabriel et al. 2004) and CIAO 4.11 (Fruscione et al. 2006). Standard filtering is applied to both MOS and pn data, i.e., up to quadruple events are retained and energy ranges of 0.2–12 keV for MOS and 0.2–15 keV for pn are considered. To remove periods of background flaring, we apply count-rate cuts of 2.1, 3, and 18 s⁻¹ and obtain effective exposure times of 58.5, 63.3, and 44.0 ks for MOS1, MOS2, and pn, respectively. We then use wavdetect on the pn data to determine the source position and its uncertainty to be (R.A., decl.[J2000]) = (09^h52^m0827 ± 003, −06°07'2780 ± 045); all errors are 1σ unless otherwise noted. With an angular resolution of 41 for pn and estimated systematic uncertainty⁵ in XMM-Newton positions of 12 (1σ), we can positively associate our X-ray point source with that of PSR J0952−0607. We generate psf images of the MOS1, MOS2, and pn data using psfgen, centered on the source position at 0.5, 1, and 3 keV, and do not find any noticeable extended emission.

Zoom In Zoom Out Reset image size

We use eregionanalyse to determine optimum regions for extracting MOS and pn source spectra. With background regions shown in Figure 1, this procedure indicated an optimal extraction radius of 10'' and background-subtracted counts of 39 ± 10, 65 ± 12, and 142 ± 20 and count rates of (0.67 ± 0.18) × 10⁻³, (1.02 ± 0.19) × 10⁻³, and (3.24 ± 0.44) × 10⁻³ s⁻¹ for MOS1, MOS2, and pn, respectively. Using epatplot, we find our spectra of PSR J0952−0607 are not affected by pile-up. We account for bad pixels and chip gaps using backscale. We then compute rmf and arf files. In order to improve statistics, we combine the MOS1 and MOS2 spectra using epicspeccombine. Spectra are binned using ftools task grppha to a minimum of 15 photons per bin for each of the combined MOS spectrum and pn spectrum.

We perform spectral fitting using Xspec 12.10.1 (Arnaud 1996). We use constant to model a possible instrumental difference between MOS and pn spectral normalizations, and we fix its value to 1 for the pn spectrum and allow it to vary for the combined MOS spectrum. To model X-ray absorption by the interstellar medium, we use tbabs with abundances from Wilms et al. (2000). To model the intrinsic spectrum of PSR J0952−0607, we consider either a single component composed of a power law (PL; powerlaw), blackbody (BB; bbodyrad), or neutron star atmosphere, or two components composed of combinations of the above. For a (nonmagnetic) neutron star hydrogen atmosphere model X-ray spectrum, we use nsatmos (Heinke et al. 2006) and fix the model parameters of neutron star mass and radius to $M=1.4\,{M}_{\mathrm{Sun}}$ and R = 10 km, respectively, and distance to d = 1.74 kpc. We also consider the nonmagnetic atmosphere model nsspec with an iron or solar composition, which are computed for fixed M = 1.4 ${M}_{\mathrm{Sun}}$ and R = 10 km and thus fixed surface gravity and gravitational redshift (Gänsicke et al. 2002).

For our first set of spectral fits, we allow the absorption parameter ${N}_{{\rm{H}}}$ to be free to vary. The results of our simultaneous fit of the pn and combined MOS spectra are given in Table 1. The top panel of Figure 2 shows the results of a spectral fit using a powerlaw model. A PL provides a generally good fit of the spectra of PSR J0952−0607, with a photon index ${\rm{\Gamma }}\approx {2.5}_{-0.4}^{+0.5}$ and unabsorbed 0.3–10 keV flux ${f}_{0.3-10}^{\mathrm{unabs}}\approx 9\times {10}^{-15}$ erg cm⁻² s⁻¹. The latter results in a luminosity L = 3 × 10³⁰ erg s⁻¹ at a distance of 1.74 kpc, which is more than 10 times lower than the upper limits determined using a short Swift exposure by Bassa et al. (2017). Even though a PL is a good fit to the data, there is possibly unmodeled excess flux at the highest energies (E ≳ 5 keV), the derived ${N}_{{\rm{H}}}\approx ({12}_{-5}^{+7})\times {10}^{20}$ cm⁻² is somewhat higher than the (4−7) × 10²⁰ cm⁻² estimated to be in the direction of PSR J0952−0607 (see Section 1), and the photon index Γ ≈ 2.5 is relatively soft and suggestive that a thermal model is possibly more appropriate. However, single component thermal models are a poor fit: BB with χ²/dof = 20/13 (see Table 1), nsatmos with χ²/dof = 18/13, and nsspec with χ²/dof = 23/13 for iron and χ²/dof = 16/13 for solar composition. The poor fits are due to the model spectra not being able to match the observed flux at high energies (≳2 keV). Note that, for the spectral fit using nsatmos, we allow the model normalization to vary, and the fit results yield values smaller than 1; while this is formally allowed and can be interpreted as emission from only a fraction of the stellar surface, it is not strictly correct because the model is computed assuming emission from the entire surface. If the normalization is fixed at one, then a much larger distance would have to be assumed (for the same inferred flux and temperature).

Zoom In Zoom Out Reset image size

Table 1.  Spectral Model Fits That Allow for a Varying ${N}_{{\rm{H}}}$

Model Fit Parameter	PL	BB	BB+PL	NSATMOS+PL	NSATMOS+PL	2BB
${N}_{{\rm{H}}}$ (1020 cm−2)	12.2−4.8+7.2	2.6−2.2+4.4	7.6−5.1+11	11−7+15	28.7−7.8+8.3	5.2−3.6+6.8
${{kT}}_{\infty }$ or ${{kT}}_{\mathrm{eff}}$ (eV)	⋯	238−33+36	175−50+44	84−41+44	38−4+2	200−33+35
${R}_{\infty }$ or ${R}_{\mathrm{em}}$ (m)a	⋯	58−17+25	110−46+180	800−500+5000	fixed 104	84−28+56
${{kT}}_{\infty }^{\mathrm{hot}}$ (eV)	⋯	⋯	⋯	⋯	⋯	1900−800+2400
${R}_{\infty }^{\mathrm{hot}}$ (m)	⋯	⋯	⋯	⋯	⋯	1.3−0.2+1.0
Γ	2.51−0.39+0.53	⋯	0.78±0.85	0.59±0.90	1.38−0.61+0.64	⋯
PL normalization (10−6)	1.74−0.35+0.46	⋯	0.34−0.25+0.37	0.26−0.18+0.56	0.76−0.39+0.56	⋯
MOS-pn normalization	1.23−0.23+0.27	1.31−0.25+0.31	1.25−0.22+0.27	1.24−0.22+0.27	1.22−0.22+0.26	1.27−0.23+0.28
${f}_{0.3-1}^{\mathrm{abs}}$ (10−15 erg cm−2 s−1)	1.9	1.9	2.0	2.0	2.1	2.0
${f}_{1-10}^{\mathrm{abs}}$ (10−15 erg cm−2 s−1)	3.6	1.3	7.6	8.0	6.3	7.1
${f}_{0.3-10}^{\mathrm{abs}}$ (10−15 erg cm−2 s−1)	5.5	3.2	9.6	10.0	8.3	9.1
χ2/dof	10.74/13	20.46/13	5.97/11	5.57/11	6.87/12	6.79/11

Notes. Xspec models: PL=powerlaw and BB=bbodyrad. For NSATMOS, $d=1.74\,\mathrm{kpc}$, M = 1.4 ${M}_{\mathrm{Sun}}$, and R = 10 km are assumed. All error bars are 1σ.

^aAssuming d = 1.74 kpc to calculate ${R}_{\infty }$ and ${R}_{\infty }^{\mathrm{hot}}$ and R = 10 km to calculate R_em.

Download table as:  ASCIITypeset image

Two component spectral models produce improved fits because a thermal component can reproduce the low energy spectrum and a PL or second hot BB reproduces the high energy emission, even above 5 keV (see Figure 2). In the cases of BB+PL and 2BB, the derived ${N}_{{\rm{H}}}$ (albeit with large uncertainty) matches that inferred from the DM of PSR J0952−0607. However, the resulting fit parameters are somewhat unusual (Γ ∼ 0.8 or ${R}_{\infty }^{\mathrm{hot}}=1\,{\rm{m}}$), and the two component model fits are not strongly preferred (f-test of BB+PL compared to PL yields a probability of 4.0% of producing by chance such a fit improvement when adding a BB). A model fit using NSATMOS+PL with the normalization of nsatmos free to vary (see above) gives results that are comparable to those of BB+PL. On the other hand, if we fix this model normalization to be unity, such that R_em = R = 10 km, then the surface temperature ${T}_{\mathrm{eff}}={4.4}_{-1.0}^{+0.5}\times {10}^{5}\,{\rm{K}}$ (at 90% confidence), and PL index Γ ≈ 1.4 is similar to that seen in other pulsars (see, e.g., De Luca et al. 2005; Guillot et al. 2016, where Γ = 1.5−2.1 is found for four pulsars), but ${N}_{{\rm{H}}}$ = (2.9 ± 0.8) × 10²¹ cm⁻² is much greater that expected. In this last case, a high ${N}_{{\rm{H}}}$ is needed to strongly reduce the low energy model flux, which is high because the nsatmos normalization is not allowed to be smaller and emission is from the entire neutron star surface. Similar results are obtained using nsatmos+nsatmos, where one component has a normalization fixed at one to model emission from the entire surface while the other component with a variable normalization could be that due to a hot spot. Also problematic with this model fit is that the temperature of the hot component hits the upper limit of the model ($\mathrm{log}{T}_{\mathrm{eff}}=6.5$); because the ratio of BB to atmosphere temperatures is seen to be ∼2 in spectral studies of other neutron stars, we expect this hot component to have a temperature that exceeds ${10}^{7}\,{\rm{K}}$. An f-test of NSATMOS+PL with fixed nsatmos normalization compared to PL yields a probability of 2.3%.

In summary, while several models can fit well the spectrum of PSR J0952−0607, none are entirely satisfactory because either the inferred ${N}_{{\rm{H}}}$ is too high or Γ is too low or too high. Nevertheless, to obtain further constraints on the surface temperature of PSR J0952−0607, we perform spectral fits using the NSATMOS+PL model with nsatmos normalization fixed at 1 and absorption column fixed at either ${N}_{{\rm{H}}}=4\times {10}^{20}$ or 1 × 10²¹ cm⁻², i.e., values that span the likely range of ${N}_{{\rm{H}}}$ (see Section 1). At these low ${N}_{{\rm{H}}}$ values compared to much higher values preferred when ${N}_{{\rm{H}}}$ is free to vary, the thermal component is forced to have a low temperature (because the emission region is the entire stellar surface), and the model fit is essentially that of a single component PL (see Table 2). From these spectral fits, we derive upper limits (at 90% confidence) on the surface temperature of 2.1 × 10⁵ and 3.1 × 10⁵ K for ${N}_{{\rm{H}}}=4\times {10}^{20}$ and 1 × 10²¹ cm⁻², respectively. Identical temperature limits are obtained when using nonmagnetic, fully ionized helium atmosphere model spectra (Ho & Lai 2001). In order to obtain a surface temperature measurement, as opposed to an upper limit, one needs ${N}_{{\rm{H}}}\,\gtrsim 1.2\times {10}^{21}$ cm⁻².

Table 2.  Spectral Model Fits with Fixed ${N}_{{\rm{H}}}$ and NSATMOS Normalization (=1)

Model Fit Parameter	PL	NSATMOS+PL	PL	NSATMOS+PL
fixed ${N}_{{\rm{H}}}$ (1020 cm−2)	4	4	10	10
${{kT}}_{\mathrm{eff}}$ (eV)	⋯	<18	⋯	<27
Γ	1.85±0.15	1.85±0.14	2.36±0.19	2.36−0.25+0.19
PL normalization (10−6)	1.22±0.17	1.22±0.17	1.62±0.22	1.62−0.22+0.20
MOS-pn normalization	1.27−0.24+0.30	1.27−0.23+0.29	1.23−0.23+0.27	1.23−0.22+0.27
${f}_{0.3-1}^{\mathrm{abs}}$ (10−15 erg cm−2 s−1)	1.6	1.6	1.9	1.9
${f}_{1-10}^{\mathrm{abs}}$ (10−15 erg cm−2 s−1)	5.3	5.3	3.9	3.9
${f}_{0.3-10}^{\mathrm{abs}}$ (10−15 erg cm−2 s−1)	6.9	6.9	5.8	5.8
χ2/dof	15.34/14	15.35/13	10.92/14	10.92/13

Note. Xspec model: PL=powerlaw. For NSATMOS, $d=1.74\,\mathrm{kpc}$, M = 1.4 ${M}_{\mathrm{Sun}}$, and R = 10 km are assumed. All error bars are 1σ, except for the 90% confidence level upper limit on ${{kT}}_{\mathrm{eff}}$.

Download table as:  ASCIITypeset image

With a small inferred size of the X-ray emitting region, such as that of a hot spot on the neutron star surface, there is the possibility of detectable X-ray pulsations if the viewing geometry is favorable and pulsations have a high enough amplitude. Unfortunately, because the time resolution of our full frame imaging mode observations is 2.6 s for MOS and 73.4 ms for pn, we are unable to search for variability due to the pulsar spin period of 1.41 ms. Nevertheless, we still perform an analysis to determine if there are other periodic signals in the data. We consider MOS2 and pn independently. We first apply a barycentric correction using barycen and DE405 ephemeris, then extract source and background 0.3–5 keV light curves using the same spatial regions as in the spectral analysis, and obtain corrected light curves using epiclccorr. The resulting light curves of PSR J0952−0607 with a bin size of 5000 s for MOS2 data and 3000 s for pn data do not show clear evidence of variability at the 6.42 hr orbital period; note that our pn and MOS exposure times span 2 and 2.5 orbital cycles, respectively. Smaller bin sizes result in time intervals where the count rate is zero. We also do not find variability using glvary, with a variability index of 0.

In this work, we analyze a recent XMM-Newton observation of the second fastest pulsar known and report on our detection of its X-ray counterpart. The data are sufficient for extraction of source spectra, and we find that these spectra can be fit well by a single PL model or a two component thermal plus PL model. We do not detect any significant variability of the source, although the XMM-Newton full frame imaging observations do not allow us to search for variations on the timescale of the short spin period of PSR J0952−0607.

With gamma-ray flux f_γ = 2.6 × 10⁻¹² erg cm⁻² s⁻¹, spin period P = 1.41 ms, and spin period time derivative $\dot{P}=4.6\times {10}^{-21}$ s s⁻¹ (Nieder et al. 2019), PSR J0952−0607 has ${f}_{\gamma }/{f}_{0.3-10}^{\mathrm{unabs}}\approx 300$, spin-down energy loss rate $\dot{E}=6.4\times {10}^{34}$ erg s⁻¹, and $L/\dot{E}\approx 5\times {10}^{-5}$. The gamma-ray to X-ray flux ratio is typical for black widow pulsars (see, e.g., Marelli et al. 2015; Salvetti et al. 2017). The $L/\dot{E}$ is well below ∼10⁻³ seen in canonical rotation-powered pulsars (Becker 2009) but similar to those seen in several millisecond pulsars in the Galactic field (see, e.g., Kargaltsev et al. 2012; Lee et al. 2018) and globular clusters (Forestell et al. 2014; Bhattacharya et al. 2017). The spin-down rate of PSR J0952−0607 could be sufficient to power its thermal X-ray luminosity. For example, deviations from beta equilibrium in the core (via rotochemical heating) of PSR J0952−0607 produce a luminosity ∼5 × 10³¹ erg s⁻¹, but this depends on uncertain properties of neutron superfluidity and proton superconductivity (Reisenegger 1995; Petrovich & Reisenegger 2010). Meanwhile, thermal creep of superfluid vortices dissipates energy at a rate of ∼(0.02−2) × 10³⁰ erg s⁻¹in PSR J0952−0607 (Alpar et al. 1984a, 1984b), and rotation-induced nuclear burning in the crust generates a heating rate of ∼5 × 10²⁹ erg s⁻¹ in PSR J0952−0607 (Gusakov et al. 2015).

In Section 3, we use model fits of the XMM-Newton spectrum of PSR J0952−0607 to constrain its surface temperature, which in turn can be used to place limits on the neutron star core temperature ${T}_{{\rm{c}}}$ through well-studied ${T}_{\mathrm{eff}}$–${T}_{{\rm{c}}}$ relations (see, e.g., Potekhin et al. 2015). For example, if the neutron star envelope is composed of iron, then

$$\begin{eqnarray}&&{T}_{{\rm{c}}}(\mathrm{Fe})=1.29\times {10}^{8}\,{\rm{K}}{\left[{g}_{14}^{-1}{\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{10}^{6}{\rm{K}}}\right)}^{4}\right]}^{5/11},\end{eqnarray} \tag{ 1 }$$

where g ≡ 10¹⁴g₁₄ cm s⁻² is surface gravity (Gudmundsson et al. 1982; Potekhin et al. 1997). On the other hand, for a fully accreted hydrogen envelope

$$\begin{eqnarray}&&{T}_{{\rm{c}}}({\rm{H}})=0.552\times {10}^{8}\,{\rm{K}}{\left[{g}_{14}^{-1}{\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{10}^{6}{\rm{K}}}\right)}^{4}\right]}^{7/17},\end{eqnarray} \tag{ 2 }$$

(Potekhin et al. 1997). Using Equations (1) and (2), our surface temperature measurement of ${T}_{\mathrm{eff}}\approx 4.4\times {10}^{5}\,{\rm{K}}$ (for varying ${N}_{{\rm{H}}}$ and thermal plus PL model spectrum) yields a core temperature of ${T}_{{\rm{c}}}=1.9\times {10}^{7}\,{\rm{K}}$ for an iron envelope or ${T}_{{\rm{c}}}=1.0\times {10}^{7}\,{\rm{K}}$ for a hydrogen envelope. Meanwhile, our surface temperature upper limit of ${T}_{\mathrm{eff}}\lt 3.1\times {10}^{5}\,{\rm{K}}$ (for fixed ${N}_{{\rm{H}}}=1\times {10}^{21}$ cm⁻² and simple PL model spectrum) yields core temperatures of ${T}_{{\rm{c}}}\lt 1.0\times {10}^{7}\,{\rm{K}}$ for an iron envelope or ${T}_{{\rm{c}}}\lt 0.56\times {10}^{7}\,{\rm{K}}$ for a hydrogen envelope.

As noted in Section 1, the r-mode instability is caused by growth of the oscillation mode through emission of gravitational waves and occurs on a timescale

$$\begin{eqnarray}&&{t}_{\mathrm{gw}}=47\,{\rm{s}}{\left(\displaystyle \frac{1000\mathrm{Hz}}{{\nu }_{{\rm{s}}}}\right)}^{6},\end{eqnarray} \tag{ 3 }$$

while the oscillation is damped by viscosity on a timescale ${t}_{\mathrm{visc}}$, which in the simplest model is due to electron–electron scattering (once the core temperature drops below the neutron superfluid transition temperature) with

$$\begin{eqnarray}&&{t}_{\mathrm{visc}}=2.2\times {10}^{5}\,{\rm{s}}{\left(\displaystyle \frac{{T}_{{\rm{c}}}}{{10}^{8}{\rm{K}}}\right)}^{2}\end{eqnarray} \tag{ 4 }$$

(Andersson & Kokkotas 2001; Shternin & Yakovlev 2008; Gusakov et al. 2014). Therefore the r-mode is unstable and grows when ${t}_{\mathrm{gw}}\lt {t}_{\mathrm{visc}}$ and is damped when ${t}_{\mathrm{gw}}\gt {t}_{\mathrm{visc}}$. This is illustrated in Figure 3, which displays the "r-mode instability window" based on the two primary parameters, ${\nu }_{{\rm{s}}}$ and ${T}_{{\rm{c}}}$, that determine ${t}_{\mathrm{gw}}$ and ${t}_{\mathrm{visc}};$ note that the dominant viscosity at ${T}_{{\rm{c}}}\gtrsim {10}^{10}\,{\rm{K}}$ is not that given by Equation (4) but by bulk viscosity. Also shown are (${\nu }_{{\rm{s}}}$, ${T}_{{\rm{c}}}$) for neutron stars in systems where these values or their limits can be determined (see Gusakov et al. 2014; Bhattacharya et al. 2017; Mahmoodifar & Strohmayer 2017; Rangelov et al. 2017; Schwenzer et al. 2017; Gonzalez-Caniulef et al. 2019, and references therein; see also Mahmoodifar & Strohmayer 2013; for simplicity, we plot fiducial ${T}_{{\rm{c}}}$ from Gusakov et al. 2014, but we note there is factor of $\lesssim 3$ uncertainty due to uncertain envelope composition and gravitational redshift). For PSR J0952−0607, we show two values of ${T}_{{\rm{c}}}$ derived from our spectral fits using either a varying or fixed ${N}_{{\rm{H}}}$ and assuming a fully accreted hydrogen envelope.

Zoom In Zoom Out Reset image size

Our results for PSR J0952−0607 can be used to inform our understanding of r-modes. Figure 3 shows that many neutron stars, including PSR J0952−0607 (see below), should be unstable to r-mode growth and thus are potentially strong sources of gravitational waves. This in part motivates r-mode searches by LIGO/Virgo (Abbott et al. 2017, 2019a; Meadors et al. 2017; Caride et al. 2019). However, spin and thermal evolution calculations indicate neutron stars should only spend a short time within the instability window and long time outside the window (Levin 1999; Heyl 2002). Thus one expects few sources to lie within the window at any one time compared to the number of sources outside the window, contrary to what is shown in Figure 3. This suggests that the instability window should be much smaller and our theoretical understanding of the physics that determines ${t}_{\mathrm{visc}}$ (e.g., neutrino emission, crust elasticity, superfluidity, and/or existence of hyperonic or quark matter; Andersson & Kokkotas 2001; Alford et al. 2012a, 2012b) is inadequate (Ho et al. 2011; Haskell et al. 2012; Chugunov et al. 2017). For example, one possible mechanism that could close the window at ${T}_{{\rm{c}}}\lt {10}^{9}\,{\rm{K}}$ is superfluid mutual friction, i.e., dissipation due to interactions between protons/electrons and superfluid neutron vortices (Haskell et al. 2009). Its possible effect on the r-mode instability window is illustrated in the left panel of Figure 4 (see Haskell et al. 2009, 2012, for details; see also Ho et al. 2011).

Zoom In Zoom Out Reset image size

As noted above, PSR J0952−0607 is inside the traditional r-mode instability region shown in Figure 3. If the measured rate of change of spin frequency ${\dot{\nu }}_{{\rm{s}}}$ is due (entirely) to energy loss by gravitational wave emission at a constant r-mode oscillation amplitude α, then the amplitude is

$$\begin{eqnarray}&&\alpha \approx 2\times {10}^{-8}{\left(\displaystyle \frac{1000\mathrm{Hz}}{{\nu }_{{\rm{s}}}}\right)}^{7/2}{\left(\displaystyle \frac{{\dot{\nu }}_{{\rm{s}}}}{{10}^{-15}\mathrm{Hz}{{\rm{s}}}^{-1}}\right)}^{1/2}\end{eqnarray} \tag{ 5 }$$

(Owen et al. 1998; Andersson & Kokkotas 2001) and is $\alpha =1\times {10}^{-7}$ for PSR J0952−0607. Another constraint on α is set by the balance between viscous heating (by the aforementioned electron–electron scattering) of an unstable r-mode and cooling by surface radiation (at these low temperatures, surface photon emission dominates core neutrino emission). This yields

$$\begin{eqnarray}\begin{array}{rcl}\alpha & = & 1.1\times {10}^{-9}{\left(\displaystyle \frac{{L}_{\mathrm{th}}}{{10}^{30}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1/2}\left(\displaystyle \frac{{T}_{{\rm{c}}}}{{10}^{7}\,{\rm{K}}}\right)\left(\displaystyle \frac{1000\,\mathrm{Hz}}{{\nu }_{{\rm{s}}}}\right)\\ & \sim & 3\times {10}^{-10}{\left(\displaystyle \frac{{L}_{\mathrm{th}}}{{10}^{30}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{31/34}\left(\displaystyle \frac{1000\,\mathrm{Hz}}{{\nu }_{{\rm{s}}}}\right)\end{array}\end{eqnarray} \tag{ 6 }$$

(Owen et al. 1998; Andersson & Kokkotas 2001), where L_th is thermal luminosity and the second equality makes use of Equation (2). Alternatively, if we assume that PSR J0952−0607 lies on the boundary between stability and instability, such that ${t}_{\mathrm{gw}}={t}_{\mathrm{visc}}$, then balance of heating and cooling yields

$$\begin{eqnarray}&&\alpha =1.6\times {10}^{-10}{\left(\displaystyle \frac{{L}_{\mathrm{th}}}{{10}^{30}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{1000\mathrm{Hz}}{{\nu }_{{\rm{s}}}}\right)}^{4}\end{eqnarray} \tag{ 7 }$$

(Mahmoodifar & Strohmayer 2013). Since we see from Figure 3 that PSR J0952−0607 is near the boundary, such that ${t}_{\mathrm{gw}}\sim {t}_{\mathrm{visc}}$, Equations (6) and (7) give similar values of $\alpha \sim 1\times {10}^{-9}$. These constraints on r-mode amplitude are among the strongest obtained thus far. For example, X-ray observations yield upper limits of $\alpha \approx {10}^{-8}-{10}^{-6}$ for millisecond pulsars and neutron stars in a low-mass X-ray binary (Mahmoodifar & Strohmayer 2013, 2017; Chugunov et al. 2017; Schwenzer et al. 2017) and of ∼2 × 10⁻⁹ for the 542 Hz pulsar 47 Tucanae aa (Bhattacharya et al. 2017).

A phenomenological approach to constraining the shape of the r-mode instability window with observed millisecond pulsars and low-mass X-ray binaries is suggested by Chugunov et al. (2017). They find observed systems require suppression of the instability in two extra regions: at ${T}_{{\rm{c}}}\sim {10}^{8}$ K and at $2\times {10}^{7}\,{\rm{K}}\lesssim {T}_{{\rm{c}}}\lesssim 3\times {10}^{7}\,{\rm{K}}$ (see the light shaded regions in right panel of Figure 4). The first region is required for consistency with the hottest neutron stars in a low-mass X-ray binary, while the second is needed for colder neutron stars in a low-mass X-ray binary and upper limits on the surface temperature of millisecond pulsars. Extension of the second constraint to even lower temperatures (see hatched regions in right panel of Figure 4) may be needed due to surface temperature limits of 47 Tucanae aa (Bhattacharya et al. 2017) and PSR J0952−0607 presented here, such that the r-mode instability is not active at ${T}_{{\rm{c}}}\lesssim 3\times {10}^{7}\,{\rm{K}}$, although isolated regions of instability are not excluded (see, e.g., Gusakov et al. 2014). Finally, Chugunov et al. (2017) argue that the r-mode instability can be almost unsuppressed at $7\,\times {10}^{7}\,{\rm{K}}\lesssim {T}_{{\rm{c}}}\lesssim {10}^{8}\,{\rm{K}}$ while still remaining consistent with observations; in such a case, this leads to a class of rapidly rotating nonaccreting neutron stars that are heated by the r-mode instability (see also Chugunov et al. 2014).

Gravitational waves from PSR J0952−0607 could be detectable if the r-mode amplitude $\alpha \gtrsim 1\times {10}^{-7}({h}_{0}/{h}_{\mathrm{sd}})$ or the pulsar has an ellipticity $\varepsilon \gtrsim 6.9\times {10}^{-10}({h}_{0}/{h}_{\mathrm{sd}})$, where h₀ is gravitational wave strain amplitude and ${h}_{\mathrm{sd}}$ is the "spin-down limit" amplitude obtained by assuming the pulsar's entire rotational energy loss is due to gravitational wave emission; note that h₀ and h_sd depend on frequency and gravitational wave mechanism (e.g., h_sd ≈ 8 × 10⁻²⁸ for an ellipticity in PSR J0952−0607). Gravitational wave searches of PSR J0952−0607 using LIGO 2015–2017 data are sensitive to h₀/h_sd ∼ 60 and thus are not able to physically constrain ε (Abbott et al. 2019b; Nieder et al. 2019).

Finally, the observing mode of the XMM-Newton data presented here is not able to resolve pulsations at the 1.41 ms spin period of PSR J0952−0607. Detection of pulsations using XMM-Newton timing mode (with a time resolution of 0.03 ms) could permit discrimination between spectral models and determination of the size of the X-ray emitting region on the neutron star. Modeling of the pulse profile could even lead to some constraints on the nuclear equation of state (e.g., Bogdanov 2013), such as that being done with observations of other millisecond pulsars using NICER (Gendreau et al. 2012).

The authors thank the anonymous referee for useful comments. W.C.G.H. appreciates use of computer facilities at the Kavli Institute for Particle Astrophysics and Cosmology. W.C.G.H. acknowledges support through grants 80NSSC19K0668 from NASA and ST/R00045X/1 from the Science and Technology Facilities Council in the United Kingdom. C.O.H. is supported by NSERC Discovery grant RGPIN-2016-04602 and a Discovery Accelerator Supplement. A.I.C. acknowledges support by RFBR according to the research project 18-32-20170. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.

Facility: XMM-Newton. -

Software: CIAO 4.11 (Fruscione et al. 2006), SAS 17.0.0 (Gabriel et al. 2004), Xspec 12.10.1 (Arnaud 1996).