NICER Detection of Thermal X-ray Pulsations from the Massive Millisecond Pulsars PSR J0740+6620 and PSR J1614-2230

We report the detection of X-ray pulsations from the rotation-powered millisecond-period pulsars PSR J0740+6620 and PSR J1614-2230, two of the most massive neutron stars known, using observations with the Neutron Star Interior Composition Explorer (NICER). We also analyze X-ray Multi-Mirror Mission (XMM-Newton) data for both pulsars to obtain their time-averaged fluxes and study their respective X-ray fields. PSR J0740+6620 exhibits a broad double-peaked profile with a separation of ~0.4 in phase. PSR J1614-2230, on the other hand, has a broad single-peak profile. The broad modulations with soft X-ray spectra of both pulsars are indicative of thermal radiation from one or more small regions of the stellar surface. We show the NICER detections of X-ray pulsations for both pulsars and also discuss the phase relationship to their radio pulsations. In the case of PSR J0740+6620, this paper documents the data reduction performed to obtain the pulsation detection and prepare for pulse profile modeling analysis.


INTRODUCTION
Rotation powered millisecond pulsars (MSPs) are usually old neutron stars (NSs) with spin periods in the ∼1-30 ms range. Their fast rotation is thought to result from the transfer of angular momentum of accreted matter from a close binary companion during a so-called "recycling" phase (Alpar et al. 1982;Radhakrishnan & Srinivasan 1982). MSPs can produce emission from radio to gamma-rays through particle acceleration and radiation processes in the pulsar magnetosphere. However, in the X-ray band substantial emission can originate from hot spots on the NS surface, and such emission can be an important probe of NS properties. This paper will focus on the X-ray emission from two pulsars found to be massive NSs.
Radio pulse timing observations of MSPs have resulted in precisely measured NS masses near or above 2 M for several objects to date. For PSR J1614−2230, Demorest et al. (2010) found M NS = 1.97 ± 0.04 M and Arzoumanian et al. (2018a) found M NS = 1.908 ± 0.016 M ; for PSR J0348+0432, Antoniadis et al. (2013) found M NS = 2.01 ± 0.04 M ; for PSR J0740+6620, Cromartie et al. (2020) found M NS = 2.14 +0.10 −0.09 M and Fonseca et al. (2021) found M NS = 2.08 ± 0.07 M . A potential fourth pulsar with high mass is PSR J1811−2405 for which Ng et al. (2020) find M NS = 2.0 +0.8 −0.5 M . Finally, radial velocity measurements in the optical band have identified potential high NS mass candidates as well. Romani et al. (2021) measured M NS = 2.13 ± 0.04 M for the MSP PSR J1810+1744, and Kandel & Romani (2020) report M NS = 2.28 +0.10 −0.09 M for the redback binary system PSR J2215+5135 (see also Linares et al. 2018). Indeed, Romani et al. (2021) explore the possibility that a substantial fraction of the NS population may exist with high masses near 2.0 M and above. How NSs can acquire such high masses is still an unresolved question, however. Depending on the evolutionary pathway of a given MSP, an appreciable amount of mass can be accreted from the companion star: for example, Tauris et al. (2011) find that PSR J1614−2230 could have accreted between ∼0.05 M and ∼0.45 M . On the other hand, Cognard et al. (2017) argue that the moderately massive NS in PSR J2222−0137 (M NS = 1.75 ± 0.06 M ) likely acquired its mass at birth and not via accretion.
Such massive NSs are of great interest for understanding the NS mass-radius relation and the equation of state (EOS) of cold, dense matter. High-mass NSs enable us to constrain the very uncertain EOS at the highest densities and connect with the EOS at lower densities obtained from less massive NSs. Furthermore, each theoretical model of NS structure predicts a particular maximum mass (Oppenheimer & Volkoff 1939; see also, e.g., Lattimer 2012), above which black hole formation occurs. For example, GW190814 indicates the merger of a black hole and a second compact object with mass 2.6 M (Abbott et al. 2020). If the second object is also a black hole, then the maximum mass of a (slowly-rotating) NS is between 2 M and 2.6 M , and there is either a small gap or no gap between the highest mass NS and lowest mass black hole. On the other hand, if knowledge gained from 2 M NSs indicates an EOS that can support much higher mass NSs, then objects such as the secondary of GW190814 could be a NS (see Abbott et al. 2018).
Even more stringent constraints on the dense matter EOS can be obtained if the mass measurement is combined with a radius estimate from modeling of the pulsed surface thermal X-ray radiation (see, e.g., Pavlov & Zavlin 1997;Miller & Lamb 1998;Bogdanov et al. 2007;Psaltis &Özel 2014). Obtaining NS radius measurements using this technique is one of the primary scientific goals of the Neutron star Interior Composition Explorer (NICER) mission. The first such constraints were presented in Miller et al. (2019) and Riley et al. (2019) based on a deep NICER exposure of the nearby isolated MSP PSR J0030+0451 .
NICER observations of massive MSPs have the potential to provide unique constraints on the properties of cold, catalyzed matter beyond nuclear density even when very X-ray faint, as discussed in depth by Miller (2016) for the case of PSR J1614−2230. Because more compact stars tend to have lower fractional modulation amplitudes, all else being equal, measurement of the thermal X-ray pulsation amplitude will provide an upper limit on the compactness GM/Rc 2 and thus a lower limit on the radius given the precisely-determined mass of this star. A lower limit to the radius of a high-mass NS, if that lower limit exceeds 10-11 km, would provide important input to nuclear theories. An upper limit on the stellar radius cannot be obtained in the same way. For example, moving the spot center toward a rotational pole can reduce the amplitude by an arbitrary factor regardless of how large the star is. A robust minimum radius can be inferred, using assumptions that otherwise maximize the overall amplitude (including the assumptions that the spot is point like and centered on the rotational equator). This is contingent, however, on being able to place constraints on any weak non-thermal pulsar emission that could potentially bias the derived radius constraints (Guillot et al. 2016).
PSR J0740+6620 is an MSP with P spin = 2.89 ms discovered in the course of the Green Bank Northern Celestial Cap (GBNCC) pulsar survey (Lynch et al. 2018). It is in a nearly circular 4.8 day binary orbit, possibly with an ultra-cool white dwarf companion (Beronya et al. 2019). Recently, Fonseca et al. (2021) obtained an updated mass measurement of M NS = 2.08 ± 0.07 M (at 68.3% confidence) for PSR J0740+6620, from the observed radio Shapiro delay, making it the NS with the largest precisely established mass (first reported by Cromartie et al. 2020). Fonseca et al. (2021) also reported an updated distance measurement from timing parallax and the Shklovskii effect of d = 1.14 +0.17 −0.15 kpc (68.3% confidence).
PSR J1614−2230 was discovered in a Parkes radio search targeting unidentified EGRET sources (Hessels et al. 2005;Crawford et al. 2006). It is an MSP with P spin = 3.15 ms bound to a massive white dwarf companion in a P b = 8.7 d orbit. The pulsar is particularly important for the NICER mission because it is one of the most massive NSs known, with the most recent measurement of M = 1.908 ± 0.016 M (Arzoumanian et al. 2018b). This precise measurement resulted from the detection of a strong signature of Shapiro delay in the binary (Demorest et al. 2010). PSR J1614−2230 has a parallax of 1.5±0.1 mas, so it lies at a parallax distance of d = 670 +50 −40 pc (Arzoumanian et al. 2018b). X-ray pulsations from this pulsar were first claimed at the ∼ 4σ confidence level based on XMM-Newton observations (Pancrazi et al. 2012). Initial exploratory X-ray observations of both pulsars with NICER, XMM-Newton, and Swift showed that they are very faint sources. We have undertaken a systematic survey with NICER of the sample of nearby rotation-powered MSPs to detect and characterize their pulsed X-ray radiation, with a focus on finding thermal X-ray pulsations that are desirable for constraints on the NS mass-radius relation and the dense matter EOS. The NICER observations of eight other MSPs were presented in Ray et al. (2019) and Guillot et al. (2019).
In the present paper, we describe NICER observations of the two nearby binary MSP systems hosting massive neutron stars, PSR J0740+6620 and PSR J1614−2230, and report the detection of thermal X-ray pulsations from both. The paper is organized as follows. In Section 2 we describe the X-ray observations and the data reduction procedures. Because the photon event data for PSR J0740+6620 have been made available for the NS radius inference analyses of Miller et al. (2021) and Riley et al. (2021), the data extraction procedure is more stringent for that pulsar than for PSR J1614−2230. In Section 3 we present our findings from the X-ray observations together with the known radio properties of these two pulsars. In Section 4, we discuss issues that are raised either by the observations themselves or the follow-on inference analyses. Finally, in Section 5 we summarize our principal conclusions.

OBSERVATIONS
The Neutron Star Interior Composition Explorer (NICER; see Gendreau & Arzoumanian 2017) is a NASA astrophysics explorer mission that has been operating on the International Space Station (ISS) since 2017 June. The NICER X-ray Timing Instrument (XTI) is an array of 52 active silicon drift detectors each paired with concentrator optics that are nominally sensitive to photons in the 0.2−12 keV range and have absolute timing precision of ∼100 nanoseconds rms for tagging of photon detection times (Prigozhin et al. 2016). The event data for both sources were processed using HEASOFT v6.28 and the NICER-specific package NICERDAS v6, together with the NICER calibration files v20200202. However, somewhat different criteria for extracting NICER event data were used for the two pulsars. Furthermore, the XMM-Newton data were obtained under a number of different circumstances, as described below.
2.1. NICER PSR J0740+6620 was observed by NICER for 1867.5 ks between 2018 September 21 and 2020 April 17, while PSR J1614−2230 was observed by NICER for 852.5 ks over the period spanning 2017 July 4 to 2021 March 1. Both pulsars are challenging to observe because of their faintness, although PSR J1614−2230 is somewhat brighter than PSR J0740+6620. Because both are faint, everything possible has been done to minimize backgrounds, including making observations primarily during ISS orbital night. Special care is warranted for PSR J0740+6620 because data sets are being used to obtain a NS radius constraint from the light curve modeling of the pulsed component (which is derived from NICER) relative to the total phase-averaged flux (which is derived from XMM-Newton imaging). Furthermore, there are nearby sources of comparable brightness that the NICER concentrating detector optics do not exclude and their contributions must be estimated. A summary of NICER observations of our target pulsars can be found in Tables 1 and 2. For the pulsar PSR J0740+6620, all events included in our and subsequent analyses were obtained when the pulsar was at least 80 • from the Sun. This minimized the effects of scattered solar X-rays and detector noise resulting from photons whose energies are just below the NICER energy range. We extract the event data in an ObsID by ObsID manner. We first filter out event data outside the 0.25−3.0 keV energy range. We also filter the observations keeping only those times that have a Cosmic Ray Cutoff Rigidity greater than 2.0 GeV/c and a space weather Kp index less than 5.0. We further require that the observations by NICER include the full complement of 52 detectors active during the observations. We exclude the detector with DET ID 34 since this detector tends to be noisy. This ensures that the event data passed to the light curve analysts have an internally consistent set of 51 out of 52 available detectors and thus utilize a constant effective area configuration.
We also demand that NICER is always pointed to better than 15 arcseconds of the source in order to include the data. The NICER vignetting function is known to be relatively flat-topped within 2 arcminutes of the instrument boresight so the pointing accuracy filter removes rare, large pointing errors as a possible significant source of brightness variation. In August 2020 maintenance on the individual detectors was begun and this involved cycling individual detectors off for a period of annealing. Thus, for our initial investigation of PSR J0740+6620, utilizing event data obtained after the detector maintenance was started would have created a NICER XTI effective area inconsistency across this date boundary and we choose to avoid this complication. The pulsar phase for each event is calculated with the photonphase tool of the PINT software package (Jing et al. 2020) using the timing solutions for each pulsar described in §2.3. The end result here is a set of good time intervals (GTIs) with event data and low background noise.
For PSR J1614−2230, because we are not yet at the stage of performing a full light curve analysis to obtain an estimated NS radius, we slightly relaxed our event acceptance criteria. Namely, we employ a Kp index constraint of 4.0 and a Sun angle constraint of 60 • . Furthermore, some of the NICER detectors suffered from a brief single event upset (SEU) on 8 July 2019, that caused them to develop anomalous time stamps until those detectors could be reset. We exclude the three days of NICER event data between the SEU occurrence and the reset. Because PSR J1614−2230 is somewhat brighter than PSR J0740+6620 this change in event acceptance criteria was found to be the best trade-off between integration time accumulation and background noise.
The NICER non-imaging field of view (FOV) with a ∼ 3 radius does not resolve X-ray sources near the targets and thus such sources can contribute to the effective background counts for the target source. Therefore, the true X-ray phase-averaged count rate for faint pulsars may be difficult to accurately determine. In order to partially mitigate this issue, we make use of XMM-Newton observations to better determine the X-ray fluxes of the pulsars and investigate the presence of other X-ray sources in the NICER FOV.

XMM-Newton
For PSR J0740+6620, we acquired three separate XMM-Newton data sets on 2019 October 26 (ObsID 0851181601), 2019 October 28 (ObsID 0851181401), and 2019 November 1 (ObsID 0851181501) through the Director's Discretionary Time program. The three EPIC instruments were operated in 'Full Frame' mode with the 'Thin' optical blocking filters in place. Due to the long instrument event readout times in 'Full Frame' mode (73.4 ms for EPIC-pn and 2.6 s for EPIC-MOS), the data cannot be used for a pulse timing analysis. As these observations occurred close to the telescope perigee passage, a large fraction of the exposures were affected by intense particle flaring. These are filtered out using the standard procedure, with special care to exclude times of high background flaring. For the EPIC-MOS cameras, we extracted events with PATTERN<=12, and for the EPIC-pn camera we extracted events with PATTERN<=4. For both instruments we excluded time ranges with significantly higher than average count rates in both the full 0.2-12 keV band and the high energy bands (10-12 keV for EPIC-MOS and 10-15 keV for EPIC-pn). This procedure resulted in the exposure times listed in Table 3. Note-All XMM-Newton exposures were obtained in 'Full-frame' mode. The "Raw" exposure times report the value of the "LIVE-TIME" FITS keyword. The "Filtered" exposure times report the sum of the GTIs after the filtering described in Section 2.2. Note-All XMM-Newton exposures were obtained in 'Full-frame' mode. The "Raw" exposure times report the value of the "LIVE-TIME" FITS keyword. The "Filtered" exposure times report the sum of the GTIs after the filtering described in Section 2.2.
For PSR J1614−2230, there are four separate observations of the pulsar in the XMM-Newton HEASARC archive. We utilize archival XMM-Newton observations (ObsID: 0404790101) originally reported by Pancrazi et al. (2012).

Radio pulsar timing
To search for pulsations, we require pulsar timing models that provide sufficient precision over the full NICER observation span to assign pulse phases to the NICER data without needing any additional trials over folding periods. These models are provided by ground-based pulsar timing programs at radio frequencies.
PSR J0740+6620 is regularly observed as part of the NANOGrav pulsar timing program (Alam et al. 2021) with both the Green Bank Telescope (GBT) and Canadian Hydrogen Intensity Mapping Experiment (CHIME). Fonseca et al. (2021) have performed an analysis of these data, extending the work done by Cromartie et al. (2020). Since this covers the entire span of the NICER data, we use that timing model to assign pulse phases for the present work.
For PSR J1614−2230, we constructed a timing model using a merged dataset from the NANOGrav 12.5yr release (Alam et al. 2021) and the Nançay radio telescope (NRT). The NRT data were obtained from MJD 57424 (2016 Feb 6) through 58921 (2020 Mar 13). The NRT backends and data processing procedure are described in Guillemot et al. (2016). We analyzed this combined data set using Tempo2 and fit a pulsar timing model containing astrometric, rotational, and binary parameters for the pulsar. The model included a time-variable dispersion measure (DM) using a second-order polynomial plus a solar wind model. The resulting pulse timing solution is presented in Table 5.  (7)  a The default definition of DM2 (and higher terms) in Tempo2 changed in 2020 June. We use the corrected Taylor series representation here.

PULSATION DETECTION
After standard processing and filtering of the NICER data as described above, we isolate the events from the 0.25 to 3.0 keV range and search for pulsations. Because for both pulsars we have established radio ephemerides we ignore trials factors and simply assign pulse phases as described above. For PSR J0740+6620 a simple Z 2 2 test yields a test statistic of Z 2 2 = 251.19, corresponding to a detection significance of 15.35σ for 574405 events in the 0.31 to 1.22 keV range. We take advantage of the fact that the filtered NICER count rate (about 0.6-1.0 s −1 ) is dominated by the background, with the pulsar contributing only a few percent to the total count rate. We can use the total count rate in each GTI interval as a proxy for the effective background count rate and employ the GTI sorting method described in Guillot et al. (2019). The GTI sorting method breaks the original GTI intervals into shorter segments and adds each segment into the evaluated time series, making sure that the overall pulsation significance increases with each addition. The lowest count rate segments are included first, and we add in segments that have increasing count rates (thus higher background) so long as the detection significance continues to increase. Once the background noise becomes sufficiently large that adding additional GTIs from the original distribution decreases rather than increases the detection significance, we stop the procedure. When we apply this method to PSR J0740+6620 we increase the overall pulsation detection significance to σ = 15.46 for 521004 events in the 0.31 to 1.18 keV range. The GTI selection is particularly effective for faint pulsars with count rates ∼ 0.01 c/s such as PSR J0740+6620. The event lists and GTIs are then saved in both text and FITS file format so that they can be utilized in subsequent processing, in particular the light curve modeling of Miller et al. (2021) and Riley et al. (2021) .
For PSR J1614−2230 we confirm the pulsations originally found by Pancrazi et al. (2012) using XMM-Newton timing-mode data. Searching the originally extracted 734.3 ks of event data for this pulsar in the 0.25 to 3.00 keV range detects the pulsations with significance Z 2 2 = 72.89 giving σ = 7.81 for 335614 events in the 0.59 to 2.20 keV range. When we optimize the good time intervals as above, our detection significance rises to σ = 8.31 but the total time is cut somewhat to 644.9 ks, the number of included events is likewise reduced to 267122, but the optimal energy range is the same at 0.59 to 2.20 keV. Figure 1 shows the folded NICER pulse profiles of PSRs J0740+6620 and J1614−2230 in the energy bands where the maximum pulse detection significance is found as determined by the Z 2 2 -test (de Jager et al. 1989;de Jager & Büsching 2010), together with their radio pulse profiles aligned in rotational phase according to the absolute pulse arrival times provided by the radio timing model of each pulsar.

DISCUSSION
The NICER observations presented here reveal PSR J0740+6620 and PSR J1614−2230 as pulsed millisecond-variable X-ray sources with phase-broadened thermal pulsations. Both NSs have masses near 2 M and thus are important X-ray sources to elucidate the cold matter EOS near the maximum NS mass. PSR J0740+6620 shows a double pulse light curve in the NICER energy range while PSR J1614−2230 shows one broad pulse. The relative phasing of the radio and X-rays shows that for PSR 0740+6620 the radio pulses are aligned with the thermal X-ray pulses and that the X-ray pulses are separated by roughly 0.4 in phase. In the case of PSR J1614-2230 the dominant radio pulse occurs near the end of the broad single thermal pulse although there is a very small second pulse in the radio light curve. Both pulsars are established as emitting from radio to γ-ray energies, indicating that non-thermal emission from the magnetosphere is present. We found no evidence, however, of non-thermal X-ray emission above 2.5 keV, almost certainly because such emission is too faint to be detected in these pulsars.
The goals of these NICER and XMM-Newton observations are to establish the existence of X-ray pulsations from both pulsars and to estimate the time-averaged flux of each pulsar. The XMM-Newton measurements help constrain  the total flux from the NS, enabling a better estimate of the pulsed emission fraction in the NICER data, and thus the estimate of the NS compactness. A correct pulsed NICER flux from the pulsar, if paired with an erroneously high pulse-averaged flux from XMM-Newton, would yield a decreased pulsed fraction and thus an incorrect lower NS radius. Conversely, a spuriously low XMM-Newton flux would yield a larger NS radius. Because of the high NS masses in these systems, even a solid lower limit in the radius of each pulsar would be valuable for constraining the NS EOS because the cores of such massive NSs contain higher-density matter. An upper limit on the compactness, along with the known NS masses, can yield a lower limit on the pulsar radius. These NICER data sets for PSR J0740+6620 and PSR J1614-2230 are very quiet in that the background count rates are strongly minimized. Also, for PSR J0740+6620 the effective area as a function of time should be very constant because we always utilized the events from 51 out of 52 active detectors.
The XMM-Newton observations of the field of PSR J0740+6620 allow us to estimate the expected count rates of the likely AGN in the field that is 3.57 arcminutes north and east from the pulsar. In fact, investigation of this X-ray point source shows that there are really two extragalactic sources within 3.157 arcseconds of one another [SDSS J074114.62+662235.2 and SDSS J074115.14+662234.9, from the Sloan Digital Sky Survey Data Release 9 (Ahn et al. 2012)] so the XMM-Newton observations we obtained can not easily distinguish if the X-rays are coming from one or both of these extragalactic sources. Using the SAS utility eregionanalyse we estimate a 0.5−3.0 keV pulsar count rate of ∼0.021 c/s for the XMM-Newton MOS1. For an assumed thermal black body pulsar spectrum with kT = 0.2 keV and converting the observed radio average dispersion measure into an estimate of the absorbing column n H = 4.5 × 10 20 cm −2 (He et al. 2013), the MOS1 count rate translates into a NICER count rate of 0.011 c/s in the 0.5−3.0 keV energy range. For the position of SDSS J074115.14+662234.9 eregionanalyse gives a NICER count rate of ∼3.1×10 −2 for this AGN and we assume a non-thermal spectrum with power law spectral index Γ ∼ 2.71. However, since the AGN (at 3.57 arcminutes away) is at the 21% level of the NICER vignetting curve (Takahashi 2013) as a function of angular radius from the look direction of the XTI, this implies the AGN will only add ∼ 0.7 of the pulsar count rate when NICER is pointed directly at PSR J0740+6620. Weak sources in the field will also contribute some counts because the NICER concentrator optics will cause some photons from these sources to fall sufficiently close to the center of the detectors to be included in the background.

CONCLUSIONS
The event data from NICER and XMM-Newton we described here for PSR J0740+6620 were used by Miller et al. (2021); Raaijmakers et al. (2021); Riley et al. (2021) to investigate the radius and EOS of the pulsar in this system. The event list we derive for this pulsar is available for download on Zenodo.org website 1 . For PSR J0740+6620 the NICER observations have shown X-ray pulsations with two broad thermal X-ray peaks separated in phase by about ∼0.4. In the case of PSR J1614−2230, we have confirmed the X-ray pulsations originally reported by Pancrazi et al. (2012) and find one broad thermal X-ray peak in the phase diagram. We show the fields of these pulsars in XMM-Newton images and find that both pulsars are faint and nearby X-ray sources will contribute significantly to their respective NICER concentrator count rates.