# Large High-precision X-Ray Timing of Three Millisecond Pulsars with NICER: Stability Estimates and Comparison with Radio

J.S. Deneva, P.S. Ray, A. Lommen, S.M. Ransom, S. Bogdanov, M. Kerr, K.S. Wood, Z. Arzoumanian, K. Black, J. Doty, et al.

► To cite this version:

J.S. Deneva, P.S. Ray, A. Lommen, S.M. Ransom, S. Bogdanov, et al., Large High-precision X-Ray Timing of Three Millisecond Pulsars with NICER: Stability Estimates and Comparison with Radio. Astrophys.J., 2019, 874 (2), pp.160. 10.3847/1538-4357/ab0966. hal-02097437

HAL Id: hal-02097437

https://hal.science/hal-02097437

Submitted on 12 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

# High-precision X-Ray Timing of Three Millisecond Pulsars with NICER: Stability Estimates and Comparison with Radio

J. S. Deneva¹, P. S. Ray², A. Lommen³, S. M. Ransom⁴, S. Bogdanov⁵, M. Kerr², K. S. Wood⁶, Z. Arzoumanian⁷, K. Black⁸, J. Doty⁹, K. C. Gendreau⁷, S. Guillot¹⁰,¹¹, A. Harding⁸, N. Lewandowska¹²,¹³, C. Malacaria¹⁴,¹⁵, C. B. Markwardt⁸, S. Price⁸, L. Winternitz⁸, M. T. Wolff², L. Guillemot¹⁶,¹⁷, I. Cognard¹⁶,¹⁷, P. T. Baker¹²,¹³, H. Blumer¹²,¹³, P. R. Brook¹²,¹³, H. T. Cromartie¹⁸, P. B. Demorest¹⁹, M. E. DeCesar²⁰, T. Dolch²¹, J. A. Ellis²², R. D. Ferdman²³, E. C. Ferrara⁸, E. Fonseca²⁴, N. Garver-Daniels¹²,¹³, P. A. Gentile¹²,¹³, M. L. Jones¹²,¹³, M. T. Lam¹²,¹³, D. R. Lorimer¹²,¹³, R. S. Lynch²⁵, M. A. McLaughlin¹²,¹³, C. Ng²⁶,²⁷, D. J. Nice²⁰, T. T. Pennucci²⁸, R. Spiewak²⁹, I. H. Stairs²⁶, K. Stovall¹⁹, J. K. Swiggum³⁰, S. J. Vigeland³⁰, and W. W. Zhu³¹

¹ George Mason University, Resident at the Naval Research Laboratory, Washington, DC 20375, USA

² U.S. Naval Research Laboratory, Washington, DC 20375, USA

³ Haverford College, Haverford, PA 19041, USA

⁴ National Radio Astronomy Observatory, Charlottesville, VA, USA

⁵ Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA

⁶ Praxis, Resident at the Naval Research Laboratory, Washington, DC 20375, USA

⁷ X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

⁸ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

⁹ Noqi Aerospace Ltd., Billerica, MA 01821, USA

¹⁰ IRAP, CNRS, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France

¹¹ Université de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France

¹² Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA

¹³ Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA

¹⁴ ST12 Astrophysics Branch, NASA Marshall Space Flight Center, Huntsville, AL 35812, USA

¹⁵ Universities Space Research Association, NSSTC, Huntsville, AL 35805, USA

¹⁶ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, LPC2E, CNRS-Université d'Orléans, F-45071 Orléans, France

¹⁷ Station de Radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France

¹⁸ University of Virginia, Department of Astronomy, P.O. Box 400325, Charlottesville, VA 22904, USA

¹⁹ National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA

²⁰ Department of Physics, Lafayette College, Easton, PA 18042, USA

²¹ Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA

²² Infinia ML, 202 Rigsbee Avenue, Durham, NC 27701, USA

²³ School of Chemistry, University of East Anglia, Norwich, NR4 7TJ, UK

²⁴ Department of Physics, McGill University, 3600 University Street, Montreal, QC H3A 2T8, Canada

²⁵ Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA

²⁶ Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada

²⁷ Dunlap Institute, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

²⁸ Hungarian Academy of Sciences MTA-ELTE Extragalactic Astrophysics Research Group, Institute of Physics, Eötvös Loránd University, Pázmány P. s. 1/A, Budapest 1117, Hungary

²⁹ Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia

³⁰ Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA

³¹ National Astronomical Observatories, Chinese Academy of Science, 20A Datun Road, Chaoyang District, Beijing 100012, People's Republic of China

Received 2018 November 13; revised 2019 February 7; accepted 2019 February 17; published 2019 April 3

## Abstract

The Neutron Star Interior Composition Explorer (NICER) is an X-ray astrophysics payload on the International Space Station. It enables unprecedented high-precision timing of millisecond pulsars (MSPs) without the pulse broadening and delays due to dispersion and scattering within the interstellar medium that plague radio timing. We present initial timing results from a year of data on the MSPs PSR B1937+21 and PSR J0218+4232, and nine months of data on PSR B1821-24. NICER time-of-arrival uncertainties for the three pulsars are consistent with theoretical lower bounds and simulations based on their pulse shape templates and average source and background photon count rates. To estimate timing stability, we use the $\sigma_z$ measure, which is based on the average of the cubic coefficients of polynomial fits to subsets of timing residuals. So far we are achieving timing stabilities $\sigma_z \approx 3 \times 10^{-14}$ for PSR B1937+21 and on the order of $10^{-12}$ for PSRs B1821-24 and J0218+4232. Within the span of our NICER data we do not yet see the characteristic break point in the slope of $\sigma_z$; detection of such a break would indicate that further improvement in the cumulative root-mean-square timing residual is limited by timing noise. We see this break point in our comparison radio data sets for PSR B1821-24 and PSR B1937+21 on timescales of >2 yr.

**Key words:** pulsars: general – pulsars: individual (PSR B1821-24, PSR B1937+21, PSR J0218+4232) – stars: neutron

**Supporting material:** tar.gz file
# 1. Introduction

The Neutron Star Interior Composition Explorer (NICER) is an X-ray instrument mounted on a movable arm on the outside of the International Space Station (ISS) and has been in operation since 2017 June (Gendreau & Arzoumanian 2017). NICER was specifically designed to study the X-ray emissions of neutron stars (NSs). The main motivation is high-precision timing of X-ray emitting pulsars, constraining the mass–radius relation of NSs, and studying their high-energy emission mechanisms. High-precision X-ray timing with NICER has already demonstrated how pulsars can be used for autonomous space navigation (Mitchell et al. 2015; Winternitz et al. 2018).

The detection of gravitational waves by LIGO (Abbott et al. 2016, 2017) has opened a new window for exploring astrophysical phenomena such as black hole and NS mergers. While LIGO is sensitive to the 10–100 Hz emission from the final inspiral of stellar-mass compact objects, pulsar timing arrays (PTAs) use long-term millisecond pulsar (MSP) timing to attempt to detect nanohertz gravitational waves from supermassive black hole binaries starting long before the system mergers. PTAs are omnidirectional, and their sensitivity improves as more observations are accumulated. The pulsars are always “on” so the experiment runs continuously, limited only by our ability to observe the pulsars. However, because PTAs have, to date, relied exclusively on radio timing, they must contend with the effects of the interstellar medium (ISM) on radio timing precision. The main ISM contributions to timing perturbations observed in PTA pulsars are from variable dispersion delays ($\propto \nu^{-2}$) and scattering delays ($\propto \nu^{-4}$) where $\nu$ is the frequency of radio emission. Most of the effect of dispersion can be corrected for during data processing. Scattering as well as variations in dispersion are stochastic processes that are difficult to model over long timescales (Shannon & Cordes 2017). They are also difficult to completely disentangle from intrinsic timing noise or from the timing perturbations caused by nHz gravitational waves.

X-ray timing observations of MSPs with NICER can help separate propagation effects from intrinsic timing noise in PTA observations. X-rays are effectively at an infinite electro-magnetic frequency compared to radio frequencies employed in PTA observations, and therefore NICER observations of MSPs are immune to the timing effects of dispersion and scattering. Ray et al. (2008) found that RXTE times of arrival (TOAs) of PSR B1821–24 agree with a GBT radio timing solution to within the X-ray TOA error bars. However, the RXTE observations produced only four TOAs within a year, which limits the timescale for estimating rotational instabilities. NICER provides superior precision, weekly observations, and the opportunity to perform this analysis for more MSPs.

Most MSPs are very faint in X-rays (e.g., Webb et al. 2004b) and many MSPs exhibit thermal X-ray pulsations that have substantially broader peaks compared to their radio counterparts (e.g., Zavlin 2006; Bogdanov & Grindlay 2009), so they are not well suited for precision timing. On the other hand, three MSPs are known to exhibit narrow nonthermal X-ray pulsations desirable for high-precision timing analyses: PSR B1821–24, PSR B1937+21, and PSR J0218+4232.

PSR B1821–24 (also known as PSR J1824–2452A) was the first radio MSP to be found in a globular cluster (Lyne et al. 1987). Pulsed X-ray emission from this isolated 3.05 ms rotator was first detected using ASCA (Saito et al. 1997); it features two remarkably narrow pulses per period with a high

pulsed fraction (~80%) and small duty cycle (Rutledge et al. 2004; Ray et al. 2008). These characteristics have made PSR B1821–24 the go-to pulsar for calibrating the absolute timing capabilities of X-ray observatories, including ASCA,<fnref n="32" /> Chandra,<fnref n="33" /> and RXTE (Rots et al. 1998).

PSR B1937+21 (PSR J1939+2134), the first MSP to be discovered (Backer et al. 1982), was found to be a pulsed X-ray source with ASCA (Takahashi et al. 2001). It was subsequently studied in more detail with BeppoSAX (Nicastro et al. 2004), RXTE (Guillemot et al. 2012), Chandra, XMM-Newton (Ng et al. 2014), and NuSTAR (Gotthelf & Bogdanov 2017). The pulse profiles of PSR B1821–24 and PSR B1937+21 show a prominent and narrow main X-ray pulse and a much weaker interpulse.

PSR J0218+4232 is a 2.32 ms pulsar (Navarro et al. 1995) bound to a white dwarf companion in a two-day binary orbit. The pulsed X-ray emission from PSR J0218+4232 has previously been studied with BeppoSAX (Mineo et al. 2000), Chandra (Kuiper et al. 2002), XMM-Newton (Webb et al. 2004a), and NuSTAR (Gotthelf & Bogdanov 2017), which revealed two moderately sharp pulses per period with a hard nonthermal spectrum.

In this paper, we use NICER data to characterize the intrinsic rotational stabilities of PSR B1821–24, PSR B1937+21, and PSR J0218+4232 and compare our findings with results from long-term radio observations from two PTA projects: the North American Nanohertz Observatory for Gravitational Waves (NANOGrav<fnref n="34" />) and the European Pulsar Timing Array (EPTA<fnref n="35" />). Section 2 describes the selection criteria we apply to NICER photons; Section 3 focuses on the pulsar timing procedure and differences between X-ray and radio timing; Section 4 compares the uncertainties of NICER pulse TOAs with simulations and theoretical predictions; Section 5 presents NICER and radio timing residuals<fnref n="36" /> for overlapping time spans; and Section 6 discusses our findings about the rotational stability of each pulsar based on NICER and radio data.

# 2. Data Selection

NICER's X-ray Timing Instrument (XTI) comprises 56 paired X-ray optics and detectors (52 currently active) coaligned to observe the same point on the sky. Each light path consists of a grazing-incidence X-ray “concentrator” and a silicon drift detector (for more detail, see Gendreau et al. 2016; Prigozhin et al. 2016). NICER is sensitive to X-rays in the 0.2–12 keV range and has a peak collecting area of 1900 cm² at 1.5 keV.

Energy deposition in theXTI’s silicon detectors, from photons or charged particles, produces an amplified charge signal. The charge signal is fed into two analog signal processing chains with different pulse shaping time constants: slow (465 ns peaking time) and fast (85 ns) (Prigozhin et al. 2016). The slow chain is optimized for energy measurement, while the fast chain is optimized for time measurement. A preset threshold in each chain produces an electronic trigger that causes the event pulse height and time stamp to be sampled

<footnote n="32">³² See https://heasarc.gsfc.nasa.gov/docs/asca/newsletters/gis_time_assign5.html.</footnote>

<footnote n="33">³³ See http://cxc.harvard.edu/contrib/arots/time/CXOClock.pdf.</footnote>

<footnote n="34">³⁴ http://nanograv.org</footnote>

<footnote n="35">³⁵ http://www.epta.eu.org</footnote>

<footnote n="36">³⁶ A residual in this case is the difference between the TOA and arrival time predicted by a model.</footnote>
and digitized. For events that trigger both chains, the fast-chain time stamp is reported; for events that trigger only the slow chain, the slow-chain time stamp is reported. In either case, the pulse height measured by the slow chain is reported. In practice, X-rays below ~1 keV trigger only the slow chain; higher-energy X-rays trigger both chains. Energetic particles, as well as gamma-rays produced as particles interact with the detector or surrounding shielding, can also cause both chains to trigger.

Event time stamps are referenced to NICER's on-board GPS receiver. They are affected by a bias (a fixed offset between a local clock and a reference clock) and an uncertainty; the latter includes contributions from time-stamping hardware accuracy, uncertainties in the lengths of cables within the instrument, and errors in the GPS receiver's realization of GPS time. Before launch, the time biases and uncertainties of the slow and fast chains were measured by laboratory equipment with calibrated biases. The fast-chain timing uncertainty was determined to be 70 ns; the difference between slow- and fast-chain time stamps was measured with 4 ns uncertainty. Thus, for all practical purposes the two analog chains are identical in terms of timing uncertainty performance. The slow- and fast-chain time biases were recorded separately for each NICER detector (typically ~250 ns for the fast chain, ~760 ns for the slow chain, with ~11 ns variations between detectors). These biases are corrected in standard processing of NICER data (the nicer-timecal routine within the HEAsoft package; see below). After correction, NICER calibrated event time stamp values refer to the time that an X-ray or particle entered the detector aperture.

Particle-induced events typically have very high amplitudes and often occur far from the center of the detector (as the 25 mm² active area of the physical detector is larger than the 2 mm diameter entrance aperture for X-rays, Prigozhin et al. 2016). They are rejected using a combination of criteria based on amplitude and offset from the detector center (Gendreau et al. 2016). The ratio (PI Ratio) of amplitudes detected by the two signal chains for the same event is strongly dependent on the offset from the detector center and provides an effective way to filter out particle background. The recorded integer PI value is the gain-corrected energy of the photon in units of 10 eV. PI_ratio=PI/PI_fast, where PI refers to the slow chain, and we exclude events with PI_ratio > 1.1 + 120/PI.

In this paper, we analyze data from 2017 July–2018 June using HEASoft 6.24³⁷ and NICERDAS 2018-04-13_V004. PSR J0218+4232 was too close to the Sun to be observed after 2018 March 30, and therefore its effective date range is 2017 July–2018 March. PSR B1821−24 also goes behind the Sun for three months out of the year, roughly November–February, so there is a gap in those data. We select Good Time Intervals (GTIs) using the following four criteria: the ISS is not within the South Atlantic Anomaly; NICER is in tracking mode; NICER is pointing within 0°015 of the source; and the source is at least 30° above the Earth's limb. We use only photons from within these GTIs.

NICER observations are typically hundreds to thousands of seconds long, and there are often multiple observations of the same source with exposure times in this range per day. For data catalog purposes, observations from a given UTC day are

grouped under an “ObsID” identifier and are downloaded, processed, and filtered together. Even after the filtering steps described above, there is still a considerable background left at low energies for some ObsIDs, often due to excessive counts in just a few detectors. This may be due to sunlight, either direct or reflected from the ISS structure, illuminating NICER detectors unevenly, combined with differences in light sensitivity between detectors. In addition, enhanced background and flaring are often observed when the ISS is in regions of low geomagnetic cutoff rigidity (“polar horn” regions), dependent on the current space weather conditions.

For each ObsID, we calculate the average count rate for each detector and exclude outlier detectors with count rates >3σ above the mean across detectors. We repeat this three times to obtain an average count rate per detector that is not contaminated by outliers. The excluded detectors as well as their total number for each ObsID differ between ObsIDs. For PSR J0218+4232, 27 of 152 ObsIDs had 1–3 excluded detectors; for PSR B1821−24, 15 of 146 ObsIDs had 1–2 excluded detectors; and for PSR B1937+21, 33 of 214 ObsIDs had 1–3 excluded detectors, and 2 ObsIDs had six excluded detectors. Finally, we filter out photons from short (8 s) stretches of data where the average counts per second are >2. This limit is based on visual inspection of diagnostic plots and is chosen to be permissive: it filters out photons from time periods affected by prominent background, which is typically in the polar horn regions. We prefer this approach over simply excluding all polar horn data because that discards a substantial amount of usable data.

For each pulsar, we derive the optimal photon energy range by computing average pulse profiles and finding the energy bounds that maximize the H-test (de Jager et al. 1989; de Jager & Büsching 2010), resulting in ranges of 0.8–6.2 keV for PSR J0218+4232, 1–5.5 keV for PSR B1821−24, and 1.15–5.55 keV for PSR B1937+21. The profiles are calculated using PINT<fnref n="38" /> with existing phase-connected radio timing solutions for all three pulsars. We exclude from further analysis photons outside these energy ranges. The resulting energy-optimized light curves are shown in Figures 1–3. The alignment between our X-ray and radio pulse profiles for the three pulsars is consistent with the literature and provides a sanity check for photon phase assignment (Cusumano et al. 2004; Knight et al. 2006; Johnson et al. 2013 for PSR J0218+4232, PSR B1821−21, and PSR B1937+21, respectively). The figures also show phaseograms of the three pulsars in grayscale, where gaps correspond to no data taken (due to scheduling priorities; in the case of PSR B1821−24, the longest gap corresponds to the pulsar being too close to the Sun to be observed), or data excluded by our selection criteria (e.g., because of high background on some days).

Overall, for PSR J0218+4232 we find 538,953 photons satisfy the selection criteria over 722.4 ks of clean exposure time, for an average of 0.75 counts s⁻¹. For PSR B1821−24, 396,435 photons and 447.5 ks of clean exposure remain after selection, with an average of 0.89 counts s⁻¹. For PSR B1937+21, we get 371,196 photons, 732.0 ks of clean exposure, and 0.51 counts s⁻¹ on average. Note that all of these count rates are total rates including source and background.

<footnote n="37">³⁷ https://heasarc.nasa.gov/lheasoft/</footnote>

<footnote n="38">³⁸ https://github.com/nanograv/pint</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top: the average pulse profile of PSR J0218+4232 from NICER data (0.80–6.20 keV) is shown in black. It is phase-aligned with an average 1484 MHz pulse profile from Nançay Radio Telescope data, shown in red, which has been scaled to the same height. Bottom: folded NICER photons vs. pulse phase and MJD. Gaps in the grayscale correspond to no data taken or data excluded by our selection criteria. Each grayscale row represents the same amount of calendar time but the exposure time per row varies. Rows with more exposure time contain more photons and therefore have smoother color variations. In both panels, two full rotations are shown for clarity, with 256 bins per rotation. A DM of 61.2365(6) pc cm⁻³ was used in calculating the delay between the frequency of the radio observations and the effective infinite frequency of NICER data.</figcaption></figure>

**3. The Timing Procedure**

To compare the timing precision and residuals of our three pulsars in X-rays and radio waves we need pulse TOAs in both energy bands spanning at least several months and overlapping in time. Table 1 lists the radio observing setups used for our three pulsars.

For PSR B1937+21, we use TOAs from the 11 yr NANOGrav data set.<fnref n="39" /> obtained from observations with the Arecibo and Green Bank telescopes (Arzoumanian et al. 2018), extended through 2018 January with additional TOAs calculated using the same data-reduction procedures and templates as the 11 yr data set. We also use TOAs from the Nançay decimetric Radio Telescope (NRT). Our combined radio TOA set for this pulsar covers ~14 yr (2004 October 14–2018 January 24).

GBT observations before 2011 used the Green Bank Astronomical Signal Processor (GASP; Demorest 2007).

Arecibo observations before 2012 used the almost identical Astronomical Signal Processor (ASP). Subsequent observations were recorded using two newer, also nearly identical backends at both observatories: the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; DuPlain et al. 2008; Ford et al. 2010) and the Puerto Rican Ultimate Pulsar Processing Instrument (PUPPI). Each backend digitized the baseband voltage signal from the receiver, performed Fourier transforms to convert the time-domain voltage signal into a channelized frequency spectrum, coherently dedispersed the data in all channels in real time, and recorded folded pulse profiles for all channels based on an existing timing solution. ASP and GASP used 4 MHz channels and recorded a folded pulse profile every 60 s. PUPPI and GUPPI used 1.56 MHz channels and recorded a folded profile every 10 s. The folded pulse profiles are summed in time and frequency, and pulse arrival times are extracted for four subbands by comparing the sum to a template using the Fourier-domain technique of Taylor (1992).

For PSR B1821−24, we use radio TOAs from observations with the NRT and Parkes telescopes. This combined TOA set spans ~7 yr (2011 August 31–2018 January 30). NRT observations were made with the Nançay Ultimate Pulsar Processing Instrument (NUPPI), which is a clone of GUPPI/PUPPI. The NUPPI setup used 4 MHz channels and recorded a folded pulse profile every 10 s. Parkes observations used both the 20 cm multibeam receiver and the upper frequency band of the coaxial 10 cm/50 cm receiver and were recorded with the fourth generation of the Parkes Digital Filterbank System, which measures the four Stokes parameters via polyphase transforms performed on Field-Programmable Gate Array processors (Manchester et al. 2013; Ferris & Saunders 2004). Data are recorded in 1024 channels and recorded to disk every 30 s, and observations are preceded by measurements of a pulsed noise diode for complex gain (polarization) calibration. For PSR J0218+4232, we use NRT (NUPPI) TOAs that also span ~7 yr (2011 August 28–2018 February 1). The NRT and Parkes TOAs we use in this work will be made publicly available in future EPTA and Parkes Pulsar Timing Array data releases. Our NICER data for the three pulsars cover the first year of the mission (2017 June 23–2018 June 30).

Pulses differ in their shape both between energy bands and from one pulse to another within each energy band. However, the average radio and X-ray pulse shapes of our three pulsars are very stable with time. This stability means that for each observing setup we can construct a noiseless template approximating the average pulse profile that is then used to extract TOAs from data taken with that setup. Our three pulsars have X-ray pulse profiles with a main pulse and interpulse ~180° apart in phase. TOA quality is maximized if we can construct a pulse shape template that fits multiple narrow resolved features. When we extract TOAs, we use a pulse shape template constructed by fitting two Gaussians to the average pulse profiles from our data span (the black NICER profiles from the top panels of Figures 1–3). Compared to using a single-Gaussian template, this also minimizes the chance that some TOAs will be calculated with respect to the main pulse while others will be calculated with respect to the interpulse. Table 2 shows the parameters of our noiseless templates for all three pulsars. In these templates, one of the peaks is centered at zero phase. In contrast, Figures 1–3 show the radio and X-ray folded pulse profiles aligned in absolute phase, which is

<footnote n="39">³⁹ The 11 yr NANOGrav TOAs and timing solution for PSR B1937+21 are available at http://data.nanograv.org.</footnote>
<table><caption>**Table 1**
Observing Setups for the Radio TOAs Used in This Paper</caption><thead><tr><th>Pulsar</th><th>MJD Range</th><th>Telescope</th><th>f<sub>center</sub><br>(GHz)</th><th>Backend</th><th>Bandwidth<br>(MHz)</th></tr></thead><tbody><tr><td>PSR J0218+4232</td><td>55801–58150</td><td>NRT</td><td>1.4</td><td>NUPPI</td><td>512</td></tr><tr><td rowspan="5">PSR B1821-24</td><td>55804–57956</td><td>NRT</td><td>1.4</td><td>NUPPI</td><td>512</td></tr><tr><td>55859–57578</td><td>NRT</td><td>2.1</td><td>NUPPI</td><td>512</td></tr><tr><td>56810–57774</td><td>NRT</td><td>2.5</td><td>NUPPI</td><td>512</td></tr><tr><td>57929–58148</td><td>Parkes</td><td>1.4</td><td>PDFB4</td><td>256</td></tr><tr><td>57929–58148</td><td>Parkes</td><td>3.0</td><td>PDFB4</td><td>1024</td></tr><tr><td rowspan="14">PSR B1937+21</td><td>53420–55974</td><td>Arecibo</td><td>1.4</td><td>ASP</td><td>64</td></tr><tr><td>53344–55968</td><td>Arecibo</td><td>2.1</td><td>ASP</td><td>64</td></tr><tr><td>56020–58312</td><td>Arecibo</td><td>1.4</td><td>PUPPI</td><td>603</td></tr><tr><td>56020–58312</td><td>Arecibo</td><td>2.1</td><td>PUPPI</td><td>460</td></tr><tr><td>53275–55243</td><td>GBT</td><td>0.8</td><td>GASP</td><td>64</td></tr><tr><td>53267–55390</td><td>GBT</td><td>1.4</td><td>GASP</td><td>48</td></tr><tr><td>55278–58302</td><td>GBT</td><td>0.8</td><td>GUPPI</td><td>186</td></tr><tr><td>55275–58301</td><td>GBT</td><td>1.4</td><td>GUPPI</td><td>642</td></tr><tr><td>55800–58367</td><td>NRT</td><td>1.4</td><td>NUPPI</td><td>512</td></tr><tr><td>55804–57574</td><td>NRT</td><td>2.1</td><td>NUPPI</td><td>512</td></tr><tr><td>56510–58301</td><td>NRT</td><td>2.5</td><td>NUPPI</td><td>512</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Same as Figure 1, but for PSR B1821-24 in the energy range 1.0–5.50 keV, and a DM of 119.8918(7) pc cm⁻³. The large gap at MJD ~58075–58160 corresponds to the pulsar being too close to the Sun to be observed.</figcaption></figure>

determined based on the radio timing solution of each pulsar. Hence the radio peak of each pulsar is at zero phase in the combined, phase-aligned plots.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same as Figure 1, but for PSR B1937+21 in the energy range 1.15–5.55 keV, and a DM of 71.01710(9) pc cm⁻³.</figcaption></figure>

While radio data are recorded as a regularly sampled time series, in X-rays we have sparse photons whose individual detection times follow Poisson statistics. Moreover, only a small percentage of the photons detected during a pulsar observation come from the pulsar as opposed to background.
<table><caption>**Table 2**
Two-Gaussian Noiseless Templates Used to Extract TOAs from NICER Data</caption><thead><tr><th>Pulsar</th><th>Phase</th><th>Amplitude</th><th>FWHM</th></tr></thead><tbody><tr><td rowspan="2">PSR J0218+4232</td><td>0.0</td><td>0.02252</td><td>0.11956</td></tr><tr><td>0.50063</td><td>0.03845</td><td>0.13552</td></tr><tr><td rowspan="2">PSR B1821-24</td><td>0.0</td><td>0.03171</td><td>0.02729</td></tr><tr><td>0.55037</td><td>0.02484</td><td>0.05216</td></tr><tr><td rowspan="2">PSR B1937+21</td><td>0.0</td><td>0.03810</td><td>0.02514</td></tr><tr><td>0.53059</td><td>0.00275</td><td>0.01177</td></tr></tbody></table>

**Note.** The phase and FWHM are in fractional phase units. The amplitude of each component is the fraction of the total counts accounted for by that component.

We use on-board GPS measurements, recorded every 10 s, to interpolate NICER's position and velocity with respect to the geocenter at each photon detection time, which is measured in the spacecraft-topocentric frame.

In order to extract an X-ray TOA from NICER data, we barycenter photon time stamps and assign a phase to each photon using PINT's photon phase routine together with the radio timing solution, which is phase-connected over the multiyear time span of the radio TOAs. We then construct a histogram of photon phases for each ObsID similar to the top panels of Figures 1–3, which show such histograms for our full NICER data span. We obtain one TOA per ObsID by using a two-Gaussian X-ray template and the maximum-likelihood method of Ray et al. (2011). The reference time for each TOA is the photon time stamp closest to the middle of the range of photon time stamps included in the ObsID. NICER spacecraft-topocentric TOAs are recorded and used together with the spacecraft's geocentric position and velocity at each TOA. This is analogous to treating ground-based TOAs, where the first transformation is from the observatory's reference frame to the geocenter and uses an ephemeris of the Earth's rotation. We include our NICER spacecraft-topocentric TOAs as supplementary electronic materials.

In our timing work we use the DE430 JPL ephemeris, Barycentric Dynamical time (TDB) units, and the TT(BIPM2015) clock realization.

For ObsIDs with total GTIs ≲100 s, there are not enough photons for the pulse peak to be significantly detected in the resulting histogram. These TOAs are excluded from the timing analysis.

## 4. TOA Uncertainties

In order to evaluate the quality of NICER TOAs, we compare the actual TOA uncertainties with estimates from theory and simulations. Even though our pulse shape templates consist of two Gaussians, as a first approximation it is useful to evaluate the single-Gaussian case analytically using only the Gaussian fitted to the higher peak. The simplest way to estimate TOA uncertainty is to calculate to what accuracy $\sigma_T$ we can localize the centroid of a Gaussian in the presence of background emission as well as statistical noise in the arrival times of source photons. This is given by the ratio of the pulse width to the pulse signal-to-noise ratio. Photon detections by NICER follow Poisson statistics; therefore, the pulse signal-to-noise in terms of total source and background photon counts within an observation is $N_{src}/\sqrt{N_{src} + N_{bkg}}$. If P is the pulse period, w is the standard deviation of the Gaussian in terms of pulse phase,

$T$ is the time span of photons yielding a single TOA, $\alpha$ is the average source photon count per second, and $\beta$ is the average background photon count per second, $\sigma_T$ is given by

$$ \sigma_T = (Pw) \sqrt{\frac{\alpha + \beta}{\alpha} \frac{1}{\sqrt{\alpha T}}} \quad (1) $$

(e.g., Lorimer & Kramer 2012).

However, all three pulsars have two-peak pulse profiles, which are better approximated as a sum of Gaussians with different amplitudes and standard deviations. Since the rate at which pulsars emit photons is phase-dependent, photon arrival times can be modeled as a nonhomogeneous Poisson process whose rate function is the two-Gaussian pulse template $h(\phi)$, where the pulse phase $\phi \in [0, 1)$ and $\int_0^1 h(\phi)d\phi = 1$. In this case, the lower limit on TOA uncertainty corresponds to the Cramer-Rao Lower Bound$^{40}$ (CRLB), as shown by Golshan & Sheikh (2007) and Winternitz et al. (2016). Golshan & Sheikh derive the CRLB as

$$ \sigma_T \ge \frac{P}{\sqrt{I_p T}}, \quad \text{where} \quad (2) $$

$$ I_p = \int_0^1 \frac{\left[ \alpha \frac{\partial}{\partial \phi} h(\phi) \right]^2}{\beta + \alpha h(\phi)} d\phi. \quad (3) $$

In the case of pulsars where $w \ll 1$ and the background is low, this can be simplified by considering only the main pulse peak and setting $\beta = 0$, so that $I_p = \alpha/w^2$.

In addition to these theoretical estimates, we simulated TOAs based on NICER instrument design parameters and average source and background photon count rates for each pulsar. We used inverse transform sampling (Devroye 1986), a technique for generating random numbers from a probability density function, to simulate pulsed events. We did this over a wide variety of integration times taking into account the predicted pulsed and unpulsed count rates and smooth two-Gaussian models of the X-ray pulse profiles, which we used as photon arrival probability density functions as a function of pulse phase. Unpulsed photons were simulated by drawing uniform random deviates over the full range of the pulse phase. Time-integrated TOAs were determined from the simulated events using the unbinned maximum-likelihood technique (Ray et al. 2011), with the smooth pulse profile model as the template, and TOA errors computed from the second moment of the resulting likelihood function.

Figures 4–6 show the uncertainties of actual NICER TOAs (where each TOA is based on data from one ObsID) and simulations. We also show estimates based on two simplified, one-Gaussian cases: Equations (1) and (2) evaluated assuming $\beta = 0$. Finally, the most realistic theoretical estimate of TOA uncertainties is from the numerically evaluated CRLB using the two-Gaussian pulse template for each pulsar and Equation (2). We find that our TOAs are consistent with the CRLB numerical lower limit as well as with simulations.

## 5. Radio and X-Ray Timing Residuals

Pulsar timing models include a parameter called dispersion measure (DM), which is the integrated column density of

<footnote>$^{40}$ The Cramer-Rao Lower Bound is a lower bound on the variance of an unbiased estimator of a fixed but unknown parameter.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** NICER TOAs vs. TOA uncertainty and exposure time for PSR J0218+4232 (black triangles) are shown together with the simulated TOAs (red) and theoretical estimates described in Section 4. The solid line corresponds to the numerical CRLB result based on the two-Gaussian pulse template (Equations (2) and (3)). Analytical results for a single Gaussian are shown with a dotted line for a simple estimate based on the ratio of pulse width and signal-to-noise (Equation (1)), and with a dashed line for the CRLB assuming zero background.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** NICER TOAs vs. TOA uncertainty and exposure time for PSR B1821-24 (black triangles) are shown together with the simulated TOAs (red) and theoretical estimates described in Section 4.</figcaption></figure>

ionized gas along the line of sight to the pulsar. The corresponding dispersion delay reflected in TOAs is $\propto$DM $\nu^{-2}$ (Lorimer & Kramer 2012). When TOAs span many years, as is the case for our radio data, the observed DM slowly and stochastically varies with time due to the changing line of sight to the pulsar and turbulence in the ISM. Any unmodeled variations manifest as red noise (i.e., noise whose spectral density is higher at lower frequencies than at high frequencies) in radio timing residuals on a timescale of months to years. X-ray TOAs are effectively at infinite electromagnetic frequency compared to radio TOAs, which means that they are immune to dispersive, diffractive, and refractive (e.g., Shannon & Cordes 2017) propagation effects due to ionized gas.

Along with any timing noise intrinsic to a pulsar, DM variations are a limiting factor for radio pulsar timing precision.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** NICER TOAs vs. TOA uncertainty and exposure time for PSR B1937+21 (black triangles) are shown together with the simulated TOAs (red) and theoretical estimates described in Section 4.</figcaption></figure>

DM variations can be modeled either by including time-derivatives of DM as additional parameters in the timing solution or by fitting offsets (DMX) to the best-fit DM at a reference epoch for TOA subsets that span a shorter period of time. DM derivatives tend to be highly covariant, unlike DMX offsets. We follow NANOGrav methods and opt for the DMX approach (Arzoumanian et al. 2018), as we can more easily extend the timing solution in the future as the data span grows.

For each pulsar, we begin with only radio TOAs and the corresponding best-fit timing solution that does not include DM derivatives or rotational frequency derivatives higher than first order. In the case of PSR B1937+21, this is the solution from the 11 yr NANOGrav data release,<fnref n="41" /> which includes DMX fits for 6 day intervals. Systematics in the NANOGrav and NRT TOA sets on PSR B1937+21 make it difficult to combine them in an unbiased manner for the purpose of measuring pulsar stability, and we treat them separately. Because both TOA sets are densely sampled and >6 yr long, two independent measurements of the stability of PSR B1937+21 are available.

We follow the procedure described in Arzoumanian et al. (2018) to generate new DMX intervals that cover the span of our radio TOAs for all three pulsars. We tailor the DMX interval length to the observation cadence for each pulsar as well as how often multifrequency radio observations occur. In the absence of multifrequency observations DMX offsets can be fitted using TOAs from several frequency subbands of the same observation. However, using TOAs from at least two different observing frequency bands that cover a wider range in frequency allows for a better DMX fit. The resulting DMX interval lengths are 6 days for PSR B1937+21, 15 days for PSR B1821-24, and 30 days for PSR J0218+4232. We disable fitting for all parameters except the DMX offsets and insert JUMPs<fnref n="42" /> every three days. The JUMPs effectively bracket each observing epoch, while allowing for epochs that may span MJDs or coordinated (usually, multifrequency) observations on adjacent days. Typical, noncoordinated observations are a week or more apart. While DMX parameters

<footnote n="41">⁴¹ https://data.nanograv.org</footnote>

<footnote n="42">⁴² A JUMP is a fitted parameter that accounts for phase offsets between data sets.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** PSR J0218+4232 timing residuals, from fits to the radio data only, for the full radio and X-ray data spans (top) and a zoomed-in range around the NICER data (bottom). NICER residuals are shown in black; NRT residuals are shown in green (1.0–1.7 GHz) and magenta (1.7–2.7 GHz).</figcaption></figure>

will preferentially absorb chromatic red noise caused by DM variations, JUMPs will preferentially absorb achromatic red noise caused by intrinsic rotational instabilities. This is equivalent to the approach taken by Taylor (1991) and Kaspi et al. (1994). However, the JUMPs and DMX offsets are still covariant. While the best fit will minimize the root-mean-square of the timing residuals, the resulting noise term separation cannot be verified to correspond to the real red noise contributions of DM variations versus intrinsic rotational instabilities by using only radio data.

Next we use Tempo<fnref n="43" /> to obtain a fit for the DMX offsets; disable fitting for the DMX offsets and enable it for all other parameters excluding DM; and remove the three-day JUMPs. After we refit again, it is the fitted rotational, astrometric, and binary parameters that absorb the red noise contribution previously absorbed by the JUMPs. Figures 7–10 show the effects of this on the residuals in the form of smoothly varying deviations from zero residual on a timescale of months to years. We include these final radio timing solutions as supplementary electronic materials.

In order to include NICER TOAs in our updated timing solution, we fit a JUMP to account for the phase difference

<figure><img src="image_9.png" /><figcaption>**Figure 8.** PSR B1821−24 timing residuals, from fits to the radio data only, for the full radio and X-ray data spans (top) and a zoomed-in range around the NICER data (bottom). NICER residuals are shown in black. NRT and Parkes residuals are shown in green (1.0–1.7 GHz), magenta (1.7–2.7 GHz), and cyan (>2.7 GHz).</figcaption></figure>

between the radio and X-ray pulse shown in Figures 1–3. Figures 7–10 show NICER and radio timing residuals with respect to the final radio timing solution from the fitting steps described above. The top panel of each figure shows the entire data span, and the bottom panel shows a zoomed area of overlap between NICER and radio TOAs. NICER TOAs are shown in black, and radio TOAs are colored according to the observing frequency. The slight remaining variations and inconsistencies between radio residuals at different observing frequencies, even after DMX offset fitting, reflect the difficulty in removing propagation effects, which is a limitation on radio timing precision.

## 6. Rotational Stability Estimates

When characterizing clock stability, the statistic of choice is typically the Allan variance, $\sigma_y^2$, which depends on second differences between clock frequency offset measurements:

$$ \sigma_y^2 = \left\langle \frac{1}{2} (\bar{y}_n - \bar{y}_{n-1})^2 \right\rangle, \quad (4) $$

where $\bar{y}_n$ is the average fractional clock frequency offset during the nth measurement interval of a certain length, and the angle

<footnote n="43">⁴³ http://tempo.sourceforge.net/</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** PSR B1937+21 timing residuals, from fits to the radio data only, for the full NRT and X-ray data spans (top) and a zoomed-in range around the NICER data (bottom). NICER residuals are shown in black. NRT residuals are shown in green (1.0–1.7 GHz), magenta (1.7–2.7 GHz), and cyan (>2.7 GHz).</figcaption></figure>

brackets denote an average over all intervals of the same length. The Allan variance is designed to quantify instability in clocks that operate at constant frequency. However, because of pulsars’ continuous energy loss, manifested as an observed spin period derivative, they act as clocks with linearly varying frequency at an a priori unknown drift rate. Therefore, the “clock noise” must be quantified in a different way.

When we obtain a best-fit timing solution for a pulsar over a span of at least one year, first-order deviations from the unknown actual pulse period are modeled and removed via fitting the pulse period, period derivative, and pulsar position on the sky. For multiyear high-precision MSP timing solutions, such as the radio fits used in this paper, we can also remove the effects of proper motion and parallax. We are interested in characterizing the remaining timing residual perturbations, whose lowest-order term is cubic. It is caused by intrinsic timing noise due, e.g., to rotational instabilities and, in the case of radio timing, imperfectly modeled propagation effects (Cordes 2013). Timing noise tends to have a “red” power spectrum, and a better measure for pulsar stability would be one that is more sensitive to red noise. Matsakis et al. (1997) introduce such a measure based on third-order variations in the timing residuals: $\sigma_z$.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** PSR B1937+21 timing residuals, from fits to the radio data only, for the full NANOGrav and X-ray data spans (top) and a zoomed-in range around the NICER data (bottom). NICER residuals are shown in black. NANOGrav residuals are shown in red (0.5–1.0 GHz), green (1.0–1.7 GHz), and magenta (1.7–2.7 GHz).</figcaption></figure>

For an interval of length $\tau$ starting at time $t_0$, we can fit a cubic polynomial to timing residuals in that interval,

$$X(t) = c_0 + c_1(t - t_0) + c_2(t - t_0)^2 + c_3(t - t_0)^3, \quad (5)$$

where $X(t)$ minimizes the sum of $[(x_i - X(t_i))/\sigma_i]^2$ over all TOAs $t_i$ with uncertainties $\sigma_i$ and residuals $x_i$. Then

$$\sigma_z \equiv \frac{\tau^2}{2\sqrt{5}} \langle c_3^2 \rangle^{1/2}, \quad (6)$$

where angle brackets denote the weighted average over the third-order coefficients of the best-fit polynomials of all nonoverlapping intervals of length $\tau$ within the set of timing residuals. In our analysis we follow the recipe for calculating $\sigma_z$ in Matsakis et al. (1997) and compute separate $\sigma_z$ values for radio and NICER TOAs for each of our three pulsars.

The full set of NANOGrav B1937+21 residuals requires special treatment: it is five years longer than the B1937+21 data set used by Matsakis et al. (1997), and a third-order polynomial fails to yield a good fit to the residuals for the largest $\tau$. Figure 11 shows third-, fourth-, and fifth-order polynomial fits for this case, and the corresponding $\sigma_z$ calculated from the third-order coefficient of each polynomial
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Full set of NANOGrav residuals for PSR B1937+21 exhibits significant red noise variations, and third- or fourth-order polynomials (dotted and dashed lines, respectively) give poor fits; a fifth-order polynomial (solid line) is necessary to obtain a good fit. In each case, the third-order coefficient is used to compute a nominal $\sigma_z$ according to Equation (6). The result for the fifth-order polynomial yields the open blue square point in Figure 15.</figcaption></figure>

fit. The resulting nominal $\sigma_z$ values differ substantially from one another.

For radio $\sigma_z$, we use the final timing solutions from Section 5, which include fitted DMX offsets. For X-ray $\sigma_z$, we adopt the best-fit astrometric parameters (position, proper motion, parallax) from these final radio solutions—the span of NICER TOAs is less than one year, too short to allow robust fitting of these parameters. ISM propagation parameters (DM, DMXs) are not relevant for X-ray timing.

We then use PINT to refit the rotational periods, period derivatives and, in the case of PSR J0218+4232, orbital parameters. This yields smaller NICER residuals compared to using the radio timing solution without refitting, because by refitting we remove some of the chromatic, imperfectly modeled ISM effects from these parameters that are present in the radio solution and do not apply to NICER TOAs. However, the new best-fit rotational period and orbital parameters retain some covariance with the fixed parameters from the initial radio solution, and therefore our updated NICER timing solution is not completely free from the influence of propagation effects plaguing the radio data nor from covariances with the astrometric parameters. NICER timing residuals with respect to this partially refitted timing solution for each pulsar are plotted in Figure 12. We include these final X-ray timing solutions for our three pulsars as supplementary electronic materials.

Figures 13–15 show $\sigma_z$ versus the length of nonoverlapping intervals $\tau$ over which we fit cubic polynomials to the radio and NICER timing residuals of PSR J0218+4232, PSR B1821−24, and PSR B1937+21. If residuals are affected only by white noise, $\sigma_z \propto \tau^{-1.5}$ (Matsakis et al. 1997); for reference, this case is illustrated by a dashed line in each figure drawn through the leftmost NICER point. A solid line shows the weighted fit for NICER $\sigma_z$ points, and a gray area shows the fit uncertainty. The redder the noise, the shallower the line slope. While the fit to NICER $\sigma_z$ points is redder than white noise for all three pulsars, the white noise slope is well within fit uncertainties. The gray region of uncertainty will shrink as we accumulate more NICER data.

For PSR B1937+21, radio $\sigma_z$ points in Figure 15 show a turn-up at $\gtrsim 600$ days in the NRT data set, and at $\gtrsim 1000$ days

<figure><img src="image_13.png" /><figcaption>**Figure 12.** NICER timing residuals with respect to the final timing solutions after refitting the rotational period, period derivative, and in the case of PSR J0218+4232, also binary parameters using only NICER TOAs, as described in Section 6.</figcaption></figure>

in the NANOGrav data set. This is roughly consistent with the results of Matsakis et al. (1997), where the turn-up occurs at $\gtrsim 800$ days. At the full length of the NANOGrav B1937+21 data set, $\sim 5000$ days, there is a spurious turn-down if $\sigma_z$ for that timescale is calculated from a third-order fit because a cubic does not describe well the residuals at that timescale (see Figure 11). Therefore, in addition to $\sigma_z$ from the poor cubic fit
<figure><img src="image_14.png" /></figure>

**Figure 13.** Measure of timing stability, $\sigma_z$ (Equation (6)), vs. the data span over which a third-order polynomial is fitted to the timing residuals of PSR J0218+4232 according to Section 6. Red and black points show $\sigma_z$ for NICER and radio data, respectively. The best fit to the NICER points is shown with a solid line and its $1\sigma$ uncertainty is shaded in gray. A dashed line plotted through the leftmost NICER point shows the slope for the case where timing precision is limited by white noise only. A vertical line marks a duration of two years.

<figure><img src="image_15.png" /></figure>

**Figure 14.** Same as Figure 13 but for PSR B1821-24.

we show (with an open blue square in Figure 15) $\sigma_z$ from the third-order coefficient of a fifth-order polynomial, which is the best fit to the full set of NANOGrav B1937+21 residuals.

The fits to NICER points are also projections of $\sigma_z$ for NICER TOAs of the three pulsars by the completion of the first 2 yr of the NICER mission, denoted by a vertical line in each figure. PSR J0218+4232 has a much wider X-ray pulse profile (Figure 1) than PSR B1821-24 and PSR B1937+21 (Figures 2 and 3, respectively), which results in larger TOA uncertainties for the same exposure time per TOA (compare Figure 4 with Figures 5 and 6) and higher $\sigma_z$ for the same timing baseline.

While the effects of white noise on timing precision can be mitigated by increasing the total number of observations as well as the total time span of TOAs, this is not true for red noise. Since both rotational and propagation effects causing red noise are stochastic and slowly varying, it becomes more prominent on large timescales. In addition, while on short timescales red noise may be subsumed in fitted timing

<figure><img src="image_16.png" /></figure>

**Figure 15.** Same as Figure 13 but for PSR B1937+21. Red, black, and blue points show $\sigma_z$ for NICER, NRT, and NANOGrav data, respectively. Third-order polynomials are used for all data spans $\tau$ and the corresponding $\sigma_z$ values are denoted by filled symbols. A third-order polynomial does not give a good fit to the full span of NANOGrav residuals (Figure 11), and results in a spurious turn-down in the NANOGrav $\sigma_z$ at $\tau \sim 5000$ days. For this case we also show with an open square the nominal $\sigma_z$ value based on the third-order coefficient of a fifth-order polynomial, which is necessary to obtain a good fit to the residuals.

parameters like rotational or DM derivatives, these parameters are stochastic and do not extrapolate well beyond the fitted data set. As a consequence, we expect the slope of $\sigma_z$ to become shallower and eventually level off with increasing $\tau$. Once that limit is reached, accumulating more TOAs no longer results in a lower cumulative root-mean-square timing residual over the entire TOA span. Figure 15 shows that this is the case for PSR B1937+21 in radio for $\tau > 10^3$ days.

NANOGrav uses Bayesian analysis first to detect whether a pulsar's residuals contain red noise and, if that is the case, to model its amplitude and spectral index (Arzoumanian et al. 2015, 2018) along with the values of other timing parameters in the presence of red noise. However, these efforts admit that some of the red noise is still absorbed, in unknown proportions, by the achromatic first period derivative and the chromatic DMX offsets. Since NICER pulsar residuals do not contain ISM-dependent chromatic red noise, our analysis can also help separate chromatic, ISM-induced red noise from achromatic red noise due to intrinsic rotational instabilities in existing long-term radio timing data sets.

Because the NICER TOA span is still $\lesssim 1$ yr, we do not yet expect to see any leveling off in NICER $\sigma_z$ values in Figures 13–15. For PSR B1821-24, the rightmost NICER $\sigma_z$ point in Figure 14 indicates that, within error bars, there may or may not be a turn-up in NICER points for $\tau \gtrsim 300$ days. This point derives from a single $c_3$ value when a cubic polynomial is fitted to the full span of NICER residuals. Accumulating more data will allow us to clarify how NICER $\sigma_z$ values behave at timescales of hundreds to thousands of days.

## 7. Discussion

One question we want to answer is: how much red noise in MSP timing data is attributable to ISM propagation effects? NICER observations do not suffer from the ISM propagation effects and consequent red noise that plagues radio observations. One of the goals of the NICER mission is to test whether
this leads to $\sigma_z$ leveling off earlier for radio TOAs, or whether intrinsic rotational instabilities dominate, in which case both radio and NICER $\sigma_z$ would level off at similar timescales $\tau$. The answer to this question has implications for detecting the stochastic, nanohertz gravitational wave background, a regime accessible to PTAs and complementary to the regime of gravitational wave emission explored by LIGO. Currently, efforts toward detecting nanohertz gravitational waves rely on ground-based radio observations of MSPs and have yet to produce a detection. Obtaining the best possible timing precision for dozens of MSPs is central to these efforts. MSP timing with NICER and future X-ray instruments may prove a valuable or perhaps even critical addition to radio PTAs, either by providing a way to better evaluate and mitigate the effects of red noise on radio PTA data sets, or by producing X-ray TOA sets on some PTA MSPs that may be used in conjunction with radio TOAs.

Radio pulsars are considered good candidates for inclusion in PTAs if TOAs with uncertainty $\lesssim 1$ µs can be obtained within integration times of $\lesssim 30$ minutes. Because X-ray photons are sparse, the NICER exposure needed to achieve the same TOA uncertainty is longer. In Figures 4–6, the best theoretical match to NICER TOA uncertainties versus exposure is the numerical CRLB result based on a smooth two-Gaussian pulse profile, denoted with a solid line. Extrapolating this to a TOA uncertainty of $1$ µs, we find that the necessary exposure is $\sim 150$ ks for PSR B1821–24, $\sim 50$ ks for PSR B1937+21, and on the order of a megasecond for PSR J0218+4232, due to its wide X-ray pulse profile, which makes a precise TOA harder to measure. However, while a radio TOA is produced from a single continuous observation recording radio flux density at a regular sampling time, X-ray TOAs may be computed from data spanning days or even weeks of observations, because the arrival time and therefore the phase of each photon is computed individually.

NICER is a technology pathfinder to future missions such as the proposed Spectroscopic Time-Resolving Observatory for Broadband Energy X-rays (STROBE-X)44 and the enhanced X-ray Timing and Polarimetry Mission (eXTP).45 We consider the implications of Equation (1) on future missions. The uncertainty of X-ray TOAs depends on the detector area in addition to integration time and background. For a given source, and for a fixed detector and mirror design, the quantity $(\alpha + \beta)/\alpha$ in Equation (1) will be the same no matter how many modules (i.e., mirror and detector) are used to make the measurement, but the rate $\alpha$ in counts/s will scale directly as area, which could be represented as $\alpha = \alpha_0(A/A_0)$, where $A_0$ is a reference area (that of NICER), $\alpha_0$ is the rate corresponding to that area, and A is the area of a hypothetical array of different size. Thus holding the integration time T constant but increasing the area will result in a more precise measurement, scaling as $1/\sqrt{A}$. However, Equation (1) could be solved for integration time T, replacing $\alpha$ and $\beta$ with rates per unit area, as

$$T = (Pw)^2 \left( \frac{\alpha_0 + \beta_0}{\alpha_0} \right) \frac{1}{\alpha_0 \sigma_T^2} \frac{A_0}{A} \quad (7)$$

Equation (7) shows that the integration time required to achieve precision of $\sigma_T$ scales as $1/A$. To get the integration time for PSR B1821–24 from $\sim 150$ ks down to 15 ks would

require increasing the array size from that of NICER ($\sim 2000$ cm$^2$) to $\sim 2$ m$^2$, maintaining the NICER design. This is less than the 5 m$^2$ area of STROBE-X. The factor by which integration time needs to be reduced depends upon the application. For detection of long-period gravitational waves it is desirable to oversample the period of the wave, while achieving the needed accuracy, $\sigma_T$, in each measurement.

It is possible to realize a substantial further benefit in X-ray arrays by increasing the ratio of mirror size to detector size. The detector radiation background is independent of the mirrors and accounts for 30%–50% of the background in the NICER design. (Other backgrounds of interest scale proportionally with mirror size.) STROBE-X and eXTP would reduce the ratio $(\alpha_0 + \beta_0)/\alpha_0$ by increasing mirror size without increasing detector size or perhaps even decreasing it, while still producing useful TOAs within a reasonable integration time. These points have been illustrated using the special case of Equation (1) but are valid more generally, as the scalings do not depend on the pulse shape.

Overall we find that NICER is performing as predicted, and we anticipate making more conclusive statements about the comparison between radio and X-ray timing stability in PSR J0218+4232, PSR B1821–24, and PSR B1937+21 when we have accumulated an additional year or more of data. Our calculations of $\sigma_z$ have demonstrated a limitation of this method for evaluating rotational stability in the case of very long sets of timing data containing red noise such as the NANOGrav B1937+21 residuals. In our future analyses of rotational stability based on NICER and radio timing residuals we anticipate using a maximum-likelihood method for estimating red noise similar to the one adopted by NANOGrav, described in Arzoumanian et al. (2015).

This work was supported by NASA through the NICER mission and the Astrophysics Explorers Program. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEA-SARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. Portions of this research performed at the Naval Research Laboratory are sponsored by NASA DPR S-15633-Y. This research has made use of the NASA Astrophysics Data System (ADS) and the arXiv. The NANOGrav Project receives support from NSF Physics Frontiers Center award number 1430284. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. The Arecibo Observatory is operated by the University of Central Florida in alliance with Ana G. Méndez-Universidad Metropolitana, and Yang Enterprises. The National Radio Astronomy Observatory and the Green Bank Observatory are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). W.W.Z. is supported by the CAS Pioneer Hundred Talents Program and the Strategic Priority Research Program of the Chinese Academy of Sciences, grant No. XDB23000000.

<footnote n="44">44 https://gammaray.nsstc.nasa.gov/Strobe-X</footnote>

<footnote n="45">45 https://www.isdc.unige.ch/expT</footnote>
Facility: NICER.

Software: astropy (Astropy Collaboration et al. 2013), PINT (Luo et al. 2019), HEAsoft (NASAHEASARC 2014), Tempo, Tempo2 (Edwards et al. 2006; Hobbs et al. 2006).

## ORCID iDs

J. S. Deneva https://orcid.org/0000-0003-1226-0793

P. S. Ray https://orcid.org/0000-0002-5297-5278

S. M. Ransom https://orcid.org/0000-0001-5799-9714

S. Bogdanov https://orcid.org/0000-0002-9870-2742

M. Kerr https://orcid.org/0000-0002-0893-4073

S. Guillot https://orcid.org/0000-0002-6449-106X

C. Malacaria https://orcid.org/0000-0002-0380-0041

L. Guillemot https://orcid.org/0000-0002-9049-8716

R. S. Lynch https://orcid.org/0000-0001-5229-7430

M. A. McLaughlin https://orcid.org/0000-0001-7697-7422

K. Stovall https://orcid.org/0000-0002-7261-594X

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, PhRvL, 116, 061102

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, PhRvL, 119, 161101

Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2015, ApJ, 813, 65

Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Backer, D. C., Kulkarni, S. R., Heiles, C., Davis, M. M., & Goss, W. M. 1982, Natur, 300, 615

Bogdanov, S., & Grindlay, J. E. 2009, ApJ, 703, 1557

Cordes, J. M. 2013, CQGra, 30, 4002

Cusumano, G., Hermsen, W., Kramer, M., et al. 2004, NuPhS, 132, 596

de Jager, O. C., & Büsching, I. 2010, A&A, 517, 9

de Jager, O. C., Raubenheimer, B. C., & Swanepoel, J. W. H. 1989, A&A, 221, 180

Demorest, P. B. 2007, PhD thesis, Univ. California

Devroye, L. 1986, Non-Uniform Random Variate Generation (New York: Springer)

DuPlain, R., Ransom, S., Demorest, P., et al. 2008, Proc. SPIE, 7019, 70191D

Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549

Ferris, R. H., & Saunders, S. J. 2004, ExA, 17, 269

Ford, J. M., Demorest, P., & Ransom, S. 2010, Proc. SPIE, 7740, 77400A

Gendreau, K., & Arzoumanian, Z. 2017, NatAs, 1, 895

Gendreau, K., Arzoumanian, Z., Adkins, P. W., et al. 2016, Proc. SPIE, 9905, 99051H

Golshan, A. R., & Sheikh, S. I. 2007, in in Proc. 63rd Annual Meeting of the Institute of Navigation (Cambridge, MA: Institute of Navigation), 413

Gotthelf, E. V., & Bogdanov, S. 2017, ApJ, 845, 159

Guillemot, L., Johnson, T. J., Venter, C., et al. 2012, ApJ, 744, 33

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

Johnson, T. J., Guillemot, L., Kerr, M., et al. 2013, ApJ, 778, 106

Kaspi, V., Taylor, J. H., & Ryba, M. F. 1994, ApJ, 428, 713

Knight, H. S., Bailes, M., Manchester, R. N., Ord, S. M., & Jacoby, B. A. 2006, ApJ, 640, 941

Kuiper, L., Hermsen, W., Verbunt, F., et al. 2002, ApJ, 577, 917

Lorimer, D. R., & Kramer, M. 2012, Handbook of Pulsar Astronomy (Cambridge: Cambridge Univ. Press)

Luo, J., Ransom, S., Demorest, P., et al. 2019, PINT: High-precision Pulsar Timing Analysis Package, Astrophysics Source Code Library, ascl:1902.007

Lyne, A. G., Brinklow, A., Middleditch, J., Kulkarni, S. R., & Backer, D. C. 1987, Natur, 328, 399

Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, 17

Matsakis, D. N., Taylor, J. H., & Eubanks, T. M. 1997, A&A, 326, 924

Mineo, T., Cusumano, G., Kuiper, L., et al. 2000, A&A, 355, 1053

Mitchell, J. W., Hassouneh, M., Winternitz, L., et al. 2015, in AIAA Guidance, Navigation, and Control Conf. (Reston, VA: AIAA), 2015-0865

Nasa High Energy Astrophysics Science Archive Research Center (HEASARC) 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Navarro, J., de Bruyn, A. G., Frail, D. A., Kulkarni, S. R., & Lyne, A. G. 1995, ApJ, 455, L55

Ng, C.-Y., Takata, J., Leung, G. C. K., Cheng, K. S., & Philippopoulos, P. 2014, ApJ, 787, 167

Nicalastro, L., Cusumano, G., Löhmer, O., et al. 2004, A&A, 413, 1065

Prigozhin, G., Gendreau, K., Doty, J. P., et al. 2016, Proc. SPIE, 9905, 99051I

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

Ray, P. S., Wolff, M. T., Demorest, P., et al. 2008, in AIP Conf. Proc. 983, 40 years of Pulsars: Millisecond Pulsars, Magnetars and More, ed. C. Bassa et al. (Melville, NY: AIP), 157

Rots, A. H., Jahoda, L., Macomb, D. J., et al. 1998, ApJ, 501, 749

Rutledge, R. E., Fox, D. W., Kulkarni, S. R., et al. 2004, ApJ, 613, 522

Saito, Y., Kawai, N., Kamae, T., et al. 1997, ApJ, 477, L37

Shannon, R. M., & Cordes, J. M. 2017, MNRAS, 464, 2075

Takahashi, M., Shibata, S., Torii, K., et al. 2001, ApJ, 554, 316

Taylor, J. H. 1991, IEEE, 79, 1054

Taylor, J. H. 1992, RSPTA, 341, 117

Webb, N. A., Olive, J.-F., & Barret, D. 2004a, A&A, 417, 181

Webb, N. A., Olive, J.-F., Barret, D., et al. 2004b, A&A, 419, 269

Winternitz, L. M. B., Mitchell, J. W., Hassouneh, M. A., et al. 2016, in 2016 IEEE Aerospace Conf. (New York: IEEE), 1

Winternitz, L. M. B., Mitchell, J. W., Hassouneh, M. A., et al. 2018, in Proc. 15th Int. Conf. on Space Operations (Reston, VA: AIAA)

Zavlin, V. E. 2006, ApJ, 638, 951