High-Precision X-ray Timing of Three Millisecond Pulsars with NICER: Stability Estimates and Comparison with Radio

TL;DR

This study demonstrates high-precision X-ray timing of millisecond pulsars using NICER, achieving stability comparable to theoretical limits and providing insights into pulsar timing noise without the interstellar medium effects present in radio observations.

Contribution

First high-precision X-ray timing analysis of millisecond pulsars with NICER, showing stability estimates and comparison with radio data to understand timing noise.

Findings

NICER achieves timing uncertainties near theoretical limits.

Timing stability of ^{-14} for PSR B1937+21.

No observed break point in _z within NICER data span.

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 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 millisecond pulsars 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 \textit{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 (RMS) 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 time scales of $> 2$ years.

Tables

Table 1. Table 1: Observing setups for the radio TOAs used in this paper.

Pulsar	MJD range	Telescope	$f_{\mathrm{center}}$	Backend	Bandwidth
			(GHz)		(MHz)
PSR J0218+4232	55801–58150	NRT	1.4	NUPPI	512
PSR B1821-24	55804–57956	NRT	1.4	NUPPI	512
	55859–57578	NRT	2.1	NUPPI	512
	56810–57774	NRT	2.5	NUPPI	512
	57929–58148	Parkes	1.4	PDFB4	256
	57929–58148	Parkes	3.0	PDFB4	1024
PSR B1937+21	53420–55974	Arecibo	1.4	ASP	64
	53344–55968	Arecibo	2.1	ASP	64
	56020–58312	Arecibo	1.4	PUPPI	603
	56020–58312	Arecibo	2.1	PUPPI	460
	53275–55243	GBT	0.8	GASP	64
	53267–55390	GBT	1.4	GASP	48
	55278–58302	GBT	0.8	GUPPI	186
	55275–58301	GBT	1.4	GUPPI	642
	55800–58367	NRT	1.4	NUPPI	512
	55804–57574	NRT	2.1	NUPPI	512
	56510–58301	NRT	2.5	NUPPI	512

Table 2. Table 2: Two-Gaussian noiseless templates used to extract TOAs from NICER data. 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.

Pulsar	Phase	Amplitude	FWHM
PSR J0218+4232	0.0	0.02252	0.11956
	0.50063	0.03845	0.13552
PSR B1821$-$24	0.0	0.03171	0.02729
	0.55037	0.02484	0.05216
PSR B1937+21	0.0	0.03810	0.02514
	0.53059	0.00275	0.01177

Equations

$${\sigma }_{\mathrm{T}}=\left(Pw\right)\sqrt{\frac{\alpha +\beta }{\alpha }}\frac{1}{\sqrt{\alpha T}}$$

$${\sigma }_{\mathrm{T}}\ge \frac{P}{\sqrt{I_{p}T}},\mathrm{where}$$

$$I_p={\int }_0^1\frac{{\left[\alpha \frac{\partial }{\partial \varphi }h\left(\varphi \right)\right]}^2}{\beta +\alpha h\left(\varphi \right)}\mathrm{d}\varphi .$$

$${\sigma }_y^2=⟨\frac{1}{2}{\left({\bar{y}}_n-{\bar{y}}_{n-1}\right)}^2⟩,$$

$$X(t)=c_0+c_1(t-t_0)+c_2(t-t_0{)}^2+c_3(t-t_0{)}^3,$$

$${\sigma }_z\equiv \frac{{\tau }^2}{2\sqrt{5}}{⟨c_3^2⟩}^{1/2},$$

$$T={\left(Pw\right)}^2\left(\frac{{\alpha }_0+{\beta }_0}{{\alpha }_0}\right)\frac{1}{{\alpha }_0{\sigma }_T^2}\frac{A_0}{A}.$$

Full text

pulsars: general — pulsars: individual (PSR B1821$-$24, PSR B1937$+$21, PSR J0218$+$4232) — stars: neutron — X-rays: stars

††software: astropy (https://ascl.net/1304.002), PINT (https://github.com/nanograv/pint), HEAsoft (https://ascl.net/1408.004), Tempo (https://ascl.net/1509.002), Tempo2 (https://ascl.net/1210.015)

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 June 2017 (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. 2018, 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 neutron star 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 time scales (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 electromagnetic 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 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 time scale 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 non-thermal 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 ($\sim$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 ASCA111See https://heasarc.gsfc.nasa.gov/docs/asca/newsletters/gis_time_assign5.html, Chandra222See http://cxc.harvard.edu/contrib/arots/time/CXOClock.pdf, 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 non-thermal 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 (NANOGrav333http://nanograv.org) and the European Pulsar Timing Array (EPTA444http://www.epta.eu.org). 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 times-of-arrival (TOAs) with simulations and theoretical predictions; Section 5 presents NICER and radio timing residuals555A residual in this case is the difference between the TOA and arrival time predicted by a model. 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) co-aligned 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 cm2 at 1.5 keV.

Energy deposition in the XTI’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 timestamp to be sampled and digitized. For events that trigger both chains, the fast-chain timestamp is reported; for events that trigger only the slow chain, the slow-chain timestamp is reported. In either case, the pulse height measured by the slow chain is reported. In practice, X-rays below $\sim 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 timestamps 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 $\sim 250$ ns for the fast chain, $\sim 760$ ns for the slow chain, with $\sim 11$ ns variations between detectors). These biases are corrected in standard processing of NICER data (the nicertimecal routine within the HEAsoft package; see below). After correction, NICER calibrated event timestamp 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 mm2 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 ${\tt PI\_RATIO}>1.1+120/{\tt PI}$.

In this paper, we analyze data from 2017 July – 2018 June using HEASoft 6.24666https://heasarc.nasa.gov/lheasoft/ and NICERDAS 2018-04-13_V004. PSR J0218+4232 was too close to the Sun to be observed after 30 March 2018, 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 (SAA); NICER is in tracking mode; NICER is pointing within $0.015^{\circ}$ 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\sigma$ 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 two 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 since 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 PINT777https://github.com/nanograv/pint 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 lightcurves are shown in Figures 1, 2, and 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 (Knight et al. 2006, Johnston et al. 2013, Cusumano et al. 2004 for PSR J0218+4232, PSR B1821$-$21, and PSR B1937+21, respectively). The figures also show phaseograms of the three pulsars in greyscale, 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*-1*. 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*-1*. For PSR B1937+21, we get 371,196 photons, 732.0 ks of clean exposure, and 0.51 counts s*-1* on average. Note that all of these count rates are total rates including source and background.

3 The Timing Procedure

To compare the timing precision and residuals of our three pulsars in X-rays and radio waves we need pulse times of arrival (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-year NANOGrav data set888 The 11-year NANOGrav TOAs and timing solution for PSR B1937+21 are available at http://data.nanograv.org 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-year data set. We also use TOAs from the Nançay decimetric Radio Telescope (NRT). Our combined radio TOA set for this pulsar covers $\sim$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 $\sim$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 multi-beam receiver and the upper frequency band of the co-axial 10cm/50cm receiver and were recorded with the fourth generation of the Parkes Digital Filterbank System (PDFB4), which measures the four Stokes parameters via polyphase transforms performed on Field-Programmable Gate Array (FPGA) processors (Manchester et al. 2013, Ferris & Saunders 2004). Data are recorded in 1024 channels and recorded to disk every 30s, 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 $\sim$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 European Pulsar Timing Array (EPTA) and Parkes Pulsar Timing Array (PPTA) 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 $\sim$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 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.

While radio data are recorded as 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. 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 photonphase routine together with the radio timing solution, which is phase-connected over the multi-year 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 $\lesssim 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_{\rm 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_{\rm src}/\sqrt{N_{\rm src}+N_{\rm 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_{\rm T}$ is given by

[TABLE]

(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 non-homogeneous Poisson process whose rate function is the two-Gaussian pulse template $h(\phi)$, where the pulse phase $\phi\in\left[0,1\right)$ and $\int_{0}^{1}h\left(\phi\right)d\phi=1$. In this case, the lower limit on TOA uncertainty corresponds to the Cramer-Rao Lower Bound999The Cramer-Rao Lower Bound is a lower bound on the variance of an unbiased estimator of a fixed but unknown parameter. (CRLB), as shown by Golshan & Sheikh (2007) and Winternitz et al. (2016). Golshan & Sheikh derive the CRLB as

[TABLE]

[TABLE]

In the case of pulsars where $w<<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: Equation 1, and Equation 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 ionized gas along the line of sight to the pulsar. The corresponding dispersion delay reflected in TOAs is $\propto{\rm 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 unmodelled 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 time scale 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. 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-year NANOGrav data release101010https://data.nanograv.org, 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 multi-frequency radio observations occur. In the absence of multi-frequency 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 JUMPs111111A JUMP is a fitted parameter that accounts for phase offsets between data sets. every three days. The JUMPs effectively bracket each observing epoch, while allowing for epochs that may span MJDs or coordinated (usually, multi-frequency) observations on adjacent days. Typical, non-coordinated observations are a week or more apart. While DMX parameters 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 vs. intrinsic rotational instabilities by using only radio data.

Next we use Tempo121212http://tempo.sourceforge.net/ 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 time scale 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 between the radio and X-ray pulse shown in Figures 1, 2, and 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^{2}_{y}$, which depends on second differences between clock frequency offset measurements:

[TABLE]

where $\bar{y}_{n}$ is the average fractional clock frequency offset during the $n$-th measurement interval of a certain length, and the angle 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 multi-year 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}$.

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

[TABLE]

where $X(t)$ minimizes the sum of $\left[(x_{i}-X(t_{i}))/\sigma_{i}\right]^{2}$ over all TOAs $t_{i}$ with uncertainties $\sigma_{i}$ and residuals $x_{i}$. Then

[TABLE]

where angle brackets denote the weighted average over the third-order coefficients of the best-fit polynomials of all non-overlapping 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 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}$ vs. the length of non-overlapping 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 grey 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 grey 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 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 time scale is calculated from a third-order fit because a cubic does not describe well the residuals at that time scale (see Figure 11). Therefore, in addition to $\sigma_{z}$ from the poor cubic fit 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 years 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 time scales. In addition, while on short time scales red noise may be subsumed in fitted timing 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 RMS 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, Arzoumanian et al. 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$ year, 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 time scales 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 time scales $\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$ $\mu$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 vs. 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 $\mu$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, since 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)131313https://gammaray.nsstc.nasa.gov/Strobe-X and the enhanced X-ray Timing and Polarimetry Mission (eXTP)141414https://www.isdc.unige.ch/extp. 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 $\left(\alpha+\beta\right)/\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}\left(A/A_{0}\right)$, 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

[TABLE]

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 cm2) to $\sim$2 m2, maintaining the NICER design. This is less than the 5 m2 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 over-sample 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 (HEASARC), 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). WWZ is supported by the CAS Pioneer Hundred Talents Program and the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000. NICER