# Timing the Pulsations of the Accreting Millisecond Pulsar SAX J1808.4–3658 during Its 2019 Outburst

Peter Bult, Deepto Chakrabarty, Zaven Arzoumanian, Keith C. Gendreau, Sebastien Guillot, Christian Malacaria, Paul S. Ray, Tod E. Strohmayer

► To cite this version:

Peter Bult, Deepto Chakrabarty, Zaven Arzoumanian, Keith C. Gendreau, Sebastien Guillot, et al.  
Timing the Pulsations of the Accreting Millisecond Pulsar SAX J1808.4–3658 during Its 2019 Out-  
burst. The Astrophysical Journal, 2020, 898 (1), pp.38. 10.3847/1538-4357/ab9827. hal-02371413

HAL Id: hal-02371413

https://hal.science/hal-02371413v1

Submitted on 6 Aug 2025

**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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Timing the Pulsations of the Accreting Millisecond Pulsar SAX J1808.4–3658 during Its 2019 Outburst

Peter Bult<sup>1,2</sup>, Deepto Chakrabarty<sup>3</sup>, Zaven Arzoumanian<sup>2</sup>, Keith C. Gendreau<sup>2</sup>, Sebastien Guillot<sup>4,5</sup>, Christian Malacaria<sup>6,7,10</sup>

<sup>1</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA

<sup>2</sup>Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>3</sup>MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

<sup>4</sup>CNRS, IRAP, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France

<sup>5</sup>Université de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France

<sup>6</sup>NASA Marshall Space Flight Center, NSSTC, 320 Sparkman Drive, Huntsville, AL 35805, USA

<sup>7</sup>Universities Space Research Association, NSSTC, 320 Sparkman Drive, Huntsville, AL 35805, USA

<sup>8</sup>Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

<sup>9</sup>Astrophysics Science Division and Joint Space-Science Institute, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

Received 2019 October 7; revised 2020 May 29; accepted 2020 May 30; published 2020 July 21

## Abstract

In this paper we present a coherent timing analysis of the 401 Hz pulsations of the accreting millisecond X-ray pulsar SAX J1808.4–3658 during its 2019 outburst. Using observations collected with the Neutron Star Interior Composition Explorer (NICER), we establish the pulsar spin frequency and orbital phase during its latest epoch. We find that the 2019 outburst shows a pronounced evolution in pulse phase over the course of the outburst. These phase shifts are found to correlate with the source flux and are interpreted in terms of hot-spot drift on the stellar surface, driven by changes in the mass accretion rate. Additionally, we find that the long-term evolution of the pulsar spin frequency shows evidence for a modulation at the Earth’s orbital period, allowing for pulsar timing based astrometry of this accreting millisecond pulsar.

*Unified Astronomy Thesaurus concepts: Millisecond pulsars (1062); Low-mass X-ray binary stars (939)*

## 1. Introduction

The low-mass X-ray binary SAX J1808.4–3658 (hereafter SAX J1808) is the canonical example of an accreting millisecond X-ray pulsar (AMXP). First discovered as an X-ray source in 1996 with the BeppoSAX satellite (in’t Zand et al. 1998), it was not until 1998 that the detection of 401 Hz pulsations with the Rossi X-ray Timing Explorer (RXTE) established this source as the first known AMXP (Wijnands & van der Klis 1998). Over the two decades since this discovery, the sample size of known AMXPs has grown substantially (Patruno & Watts 2012; Campana & di Salvo 2018). Despite this growth, however, SAX J1808 remains the cornerstone of its class. This status is in part due to its relatively close proximity, at an estimated distance of 3.5 kpc (Galloway & Cumming 2006), and in part due to its short, semiregular outburst recurrence times of 2–4 yr. Together, these qualities have enabled extensive investigations of its pulsations (Poutanen & Gierliński 2003; Hartman et al. 2008; Burderi et al. 2009; Patruno et al. 2017; Sanna et al. 2017), its aperiodic variability (Wijnands et al. 2001; Patruno et al. 2009b; Bult & van der Klis 2015a, 2015b), its spectral characteristics (Gierliński et al. 2002; Cackett et al. 2009; Papitto et al. 2009; di Salvo et al. 2019), and its thermonuclear X-ray bursts (in’t Zand et al. 2001; Chakrabarty et al. 2003; Galloway & Cumming 2006; Bhattacharyya & Strohmayer 2007; in’t Zand et al. 2013).

The characteristic property of an AMXP is its coherently pulsed emission. These oscillations have nearly sinusoidal profiles that are attributed to a localized emission region, a hot spot, on the neutron star surface. Such a hot spot arises due to

the dynamically relevant stellar magnetosphere, which channels in-falling plasma, so that the accretion column impacts the stellar surface in confined regions near the magnetic poles. These pulsations reveal a host of information about the neutron star, with the precise waveform encoding information on the stellar compactness (Poutanen & Gierliński 2003), and their phase evolution with time revealing the spin evolution of the binary and star itself (Hartman et al. 2008).

For SAX J1808 specifically, coherent timing analyses of the pulsations revealed that the X-ray outbursts do not show a measurable accretion-powered spin change (Hartman et al. 2008, 2009; Patruno et al. 2012, 2017; Sanna et al. 2017), while between outbursts the source spins down at a steady rate of ≈1.5 × 10⁻¹⁵ Hz s⁻¹. Additionally, these analyses revealed the 2.01 hr binary orbit (Chakrabarty & Morgan 1998), and established a peculiar long-term variation of the binary period. The first decade of monitoring showed that the binary period was expanding at an anomalously fast rate of ≈3.5 × 10⁻¹² s s⁻¹ (di Salvo et al. 2008; Hartman et al. 2008), which was explained as being due to highly non-conservative mass-transfer (di Salvo et al. 2008; Burderi et al. 2009), or, alternatively, due to spin-orbit angular momentum exchange in the companion star (Hartman et al. 2008; Patruno et al. 2012), akin to short-term orbital variability observed in millisecond “black widow” radio pulsars (Arzoumanian et al. 1994; Doroshenko et al. 2001). Later outbursts, in 2011 and 2015, showed that the rate of orbital expansion is itself evolving rapidly. These measurements have been analyzed in great detail (Patruno et al. 2017; Sanna et al. 2017), but for a lack of data points, no single conclusion could yet be reached.

On 2019 July 30, Russell et al. (2019b) reported a brightening in the optical i'-band from the direction of SAX J1808, which is often a precursor to a full X-ray outburst

<footnote n="10"><sup>10</sup> NASA Postdoctoral Fellow.</footnote>
(Russell et al. 2019a). Although this predictor had not previously been tested for SAX J1808 specifically, various X-ray, UV, and radio observatories were triggered to monitor SAX J1808 over the following days. On 2019 August 6 a further increase in the optical i'-band was observed, coincident with a flux increase in the X-ray band (Goodwin et al. 2019b) and UV (Parikh & Wijnands 2019), suggesting the source was indeed starting a new outburst. On 2019 August 7 the X-ray flux started increasing rapidly, and the detection of 401 Hz pulsations (Bult et al. 2019c) with the Neutron Star Interior Composition Explorer (NICER; Gendreau & Arzoumanian 2017) confirmed that SAX J1808 was indeed active again.

Following confirmation of a new outburst, we executed an extensive, high cadence monitoring campaign with NICER. In this paper we report on a timing analysis of the coherent pulsations for these data.

## 2. Observations

We observed SAX J1808 with NICER between 2019 July 30 and 2019 August 31 for a total unfiltered exposure of 355.4 ks. These data are available under program/target numbers 205026 and 258401. At the time of the final observation included in the present analysis, the source had completed the high luminosity phase of its outburst cycle, but remained visible in its characteristic prolonged flaring tail (Wijnands et al. 2001; Patruno et al. 2009b, 2016; Bult et al. 2019b).

We processed all data using NICERDAS version 6a, which is distributed as part of HEASOFT v6.26. We applied standard screening criteria, limiting our analysis to time intervals with a pointing offset <54", a bright Earth limb angle >30°, a dark Earth limb angle >15°, and collected outside the South Atlantic Anomaly. The standard processing also applies a background screening filter, which rejects all time intervals where the rate of saturating particle events (overshoots) is greater than 1 ct s⁻¹ detector⁻¹, or greater than 1.52 × COR_SAX⁻⁰.⁶³³, where COR_SAX<fnref n="11" /> gives the cutoff rigidity of the Earth's magnetic field, in units of GeV c⁻¹. For the SAX J1808 observations, we find that this approach is often too conservative, as statistical fluctuations in the overshoot rate are leading to many unnecessary 1–10 s gaps in the light curve. We therefore applied a smoothing average filter to the overshoot light curve using a 5 second window. Additionally, we relaxed the cutoff threshold for the overshoots to 1.5 ct s⁻¹ detector⁻¹ and increased the scaling factor of the COR_SAX expression from 1.52 to 2.0. A visual comparison of light curves based on standard and relaxed filtering showed that no spurious signals were introduced in this way. Using these criteria, we retained 257.7 ks of good time exposure. These good time data were corrected to the solar system barycenter using the BARYCORR tool, where we used the optical coordinates of Hartman et al. (2008) and the JPL DE405 planetary ephemeris (Standish 1998). All times reported in this work are therefore expressed in units of Barycentric Dynamical Time. Finally, we estimated the background contributions to our data from NICER observations of the RXTE blank-field regions (Jahoda et al. 2006).

Two thermonuclear X-ray bursts were observed with NICER, the first on August 9 and the second on August 21. We exclude both events from the analysis presented here. An initial analysis of the second bursts is presented in Bult et al. (2019a). An analysis of the first burst, and a separate, detailed review of the spectroscopic and stochastic variability will be presented elsewhere.

## 3. Coherent Timing

We performed a coherent timing analysis of the 401 Hz pulsations of SAX J1808. For each continuous good time interval we corrected the photon arrival times for the Doppler delays of the binary motion based on the preliminary ephemeris of Bult et al. (2019c). We then folded the data on the pulsar period to construct a pulse profile. This profile was decomposed into its Fourier components by fitting a constant plus k harmonic sinusoids. Each sinusoid thus has a fixed frequency $v_k = k v_p$, where $v_p$ is the pulsar spin frequency. Additionally, for each sinusoid we determined the phase, $\phi_k$, and amplitude, $A_k$, from which we derive a fiducial pulse arrival time and a fractional sinusoidal amplitude, r, as

$$r = \frac{A_k}{N_{\gamma} - B}, \quad (1)$$

where $N_γ$ gives the number of photons in the pulse profile, and B gives the estimated number of photons contributed by the background emission. A harmonic was deemed significant if its measured amplitude exceeded the 99% confidence amplitude threshold of the noise distribution. When a harmonic was not significantly detected, we instead computed an upper limit on the amplitude as the minimum signal that would have produced a measured amplitude in excess of the noise threshold 95% of the time.

Pulsations at the fundamental frequency are detected at high significance throughout the outburst (Figure 1, first and second panel). For a number of segments, however, the integration time was too short to reach a detection. Given the small number of cases, we did not attempt to combine these segments with neighboring intervals; we show an upper limit on the amplitude instead. The second harmonic was rarely detected, mainly because its amplitude is small below 3 keV (Patruno et al. 2009a), where NICER is most sensitive. Indeed, by accumulating more exposure in each pulse profile, we could detect the second harmonic significantly in more segments. However, the required integration times were long (>5 ks), greatly reducing the number of measured pulse arrival times. Hence, we focus our analysis on the arrival times of the fundamental pulsation.

The series of pulse arrival times was modeled as a constant spin frequency and a circular Keplerian orbit using TEMPO2 (Hobbs et al. 2006). Hence, our initial model consisted of four parameters: the binary orbital period, $P_b$, the projected semimajor axis, $a_x \sin i$, the time of passage through the ascending node, $T_{asc}$, and the spin frequency, $v_p$. After an initial fit, we repeated the epoch-folding method using the refined ephemeris, and repeated the entire procedure until the solution converged.

At a reduced $\chi^2$ (χ²) of 7.2 for 259 degrees of freedom (dof), we find that our initial model provided a poor description to the set of arrival times, leaving an approximately parabolic trend in the phase residuals (Figure 1, third panel). This trend suggests that the spin frequency is changing over the course of the

<footnote n="11">¹¹ The COR_SAX parameter is based on a model for the cutoff rigidity that was originally developed for the BeppoSAX satellite. It has no relation to the source SAX J1808.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Temporal evolution of the 2019 outburst of SAX J1808. Top panel: the 0.5–10 keV light curve. Second panel: pulse fractional sinusoidal amplitudes, with solid triangles giving the upper limits for nondetections of the fundamental component, and open triangles the upper limits for the second harmonic. Third through fifth panels: pulse phase residuals relative to the indicated phase model (see also Table 1).</figcaption></figure>

outburst. Indeed, including a spin frequency derivative in our timing model improves the fit statistic to a $\chi_r^2$ of 2.27 for 258 dof. Considering the fit residuals, however, we find large phase residuals near the start and end of our data set (Figure 1, fourth panel). This suggests that additional, higher spin frequency derivatives, for which there is little physical motivation, are needed to describe the data

An important consideration in the coherent timing analysis of accreting millisecond pulsars, however, is that the pulse profiles can change shape, leading to both gradual and abrupt shifts in the observed phase (see Patruno & Watts 2012 for a

review). Indeed, in SAX J1808 the pulse waveform is known to be sensitive to the source flux (Patruno et al. 2009c). We therefore generalized our timing model by expressing the pulsar phase as

$$ \varphi(t, F) = \varphi_0 + \nu_p t + \varphi_{\text{orb}}(t) + \varphi_{\text{bias}}(F_{\text{bol}}), \quad (2) $$

where $\varphi_0$ gives a reference phase, $\nu_p t$ gives the linear term due to the constant spin frequency, $\varphi_{\text{orb}}(t)$ gives the phase correction associated with the binary orbit, and $\varphi_{\text{bias}}(F_{\text{bol}})$ gives a phase-bias correction term as a function of bolometric flux, $F_{\text{bol}}$. In order to specify a bias function, we note that numerical simulations indicate that the azimuthal position of the hot spot is proportional to the magnetospheric radius (Kulkarni & Romanova 2013). Additionally, the magnetospheric radius is predicted to scale with the mass accretion rate as $\dot{M}^{-1/5}$ (D’Angelo & Spruit 2010; Kulkarni & Romanova 2013). Since $\dot{M}$ should be proportional to the bolometric flux, we define a simple power-law model

$$ \varphi_{\text{bias}}(F_{\text{bol}}) = b F_{\text{bol}}^{\Gamma}, \quad (3) $$

with $\Gamma$ fixed at $-1/5$, and $b$ a free scaling factor, which we will refer to as the scale of the bias. Finally, for simplicity, we use count rate as an approximate substitute for the flux. Applying this generalized model to our arrival times, we obtain a $\chi_r^2$ of 2.08 for 258 dof. While this fit statistic is still not formally acceptable, it performs as well as the quadratic phase model. Furthermore, we note that phase residuals for this model (Figure 1, fifth panel) are no longer structured, such that the poor statistic can be attributed to residual timing noise (Patruno & Watts 2012). Similar results are obtained by assuming the power-law index scales with the Alfvén radius (see Section 6.2). That is, setting $\Gamma = -2/7$ yields $b = -0.83 \pm 0.03$ for $\chi_r^2 = 2.16$. The spin frequency obtained from this fit is $6 \times 10^{-9}$ Hz smaller than that obtained from assuming an index of $1/5$. Additionally, by setting $\Gamma = 1$ we can replicate the approach of Patruno et al. (2009c). This model yields a bias of $b = (6.1 \pm 0.2) \times 10^{-4}$, with a spin frequency that is $4 \times 10^{-9}$ Hz smaller than one obtained using an index of $1/5$. At a $\chi_r^2 = 5.0$ this model does better than the linear phase model; nonetheless, it is strongly disfavored relative to the physically motivated model flux-adjusted phase model adopted here.

Finally, we note that our data quality is sufficient that the Keplerian parameters describing the binary orbit are effectively decoupled from the parameters describing the spin frequency evolution. Indeed, in all three models discussed above, the Keplerian parameters agree within their respective uncertainties. Hence, in Table 1 we list the best-fit timing solution for the orbit separately from spin frequency parameters obtained from the three phase models, with the orbital parameter values being those obtained using the linear phase model.

## 4. Pulse Energy Dependence

To investigate the photon-energy dependence of the pulsations, we divided the data into 25 energy bins, optimized such that each bin contains an approximately equal number of photons. This dynamically binned pulse energy spectrum was calculated using four representative time intervals. The first
<table><caption>**Table 1**
Timing Solution for the 2019 Outburst of SAX J1808</caption><thead><tr><th>Parameter</th><th>Value</th><th>Uncertainty</th></tr></thead><tbody><tr><td>Epoch (MJD)</td><td>58715</td><td>...</td></tr><tr><td>a<sub>x</sub> sin i (lt-ms)</td><td>62.8091</td><td>7.4 × 10<sup>-3</sup></td></tr><tr><td>P<sub>b</sub> (s)</td><td>7249.1552</td><td>3.0 × 10<sup>-3</sup></td></tr><tr><td>T<sub>asc</sub> (MJD)</td><td>58715.0221031</td><td>2.0 × 10<sup>-6</sup></td></tr><tr><td>ε</td><td>&lt;3 × 10<sup>-4</sup></td><td></td></tr><tr><td colspan="3">Linear phase model</td></tr><tr><td>ν (Hz)</td><td>400.975209741</td><td>1.2 × 10<sup>-8</sup></td></tr><tr><td>χ<sup>2</sup>/dof</td><td>1854.79/259</td><td></td></tr><tr><td colspan="3">Quadratic phase model</td></tr><tr><td>ν (Hz)</td><td>400.975210055</td><td>2.8 × 10<sup>-8</sup></td></tr><tr><td>ν̇ (Hz s<sup>-1</sup>)</td><td>-3.02 × 10<sup>-13</sup></td><td>1.3 × 10<sup>-14</sup></td></tr><tr><td>χ<sup>2</sup>/dof</td><td>585.30/258</td><td></td></tr><tr><td colspan="3">Flux-adjusted phase model</td></tr><tr><td>ν (Hz)</td><td>400.975209833</td><td>1.0 × 10<sup>-8</sup></td></tr><tr><td>b</td><td>-0.87</td><td>0.03</td></tr><tr><td>Γ</td><td>-1/5</td><td>fixed</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>535.95/258</td><td></td></tr></tbody></table>

**Note.** Uncertainties give the 1σ statistical errors and the upper limit is quoted at the 95% c.l. The $\epsilon$ parameter gives the orbital eccentricity.

interval covers August 10 through 17, which is when the outburst peaked and the broadband fractional pulse amplitude remained constant at 5%. The second interval covers August 18 through 20, which is where the source flux showed a rebrightening with respect to its long-term decay, and, likewise, the broadband fractional pulse amplitude temporarily increased to 7%. The third interval covers August 21 through 25, where the outburst continued its decay and the fractional pulse amplitude varied between 4% and 7%. Finally, the fourth group covers August 26 through 31 where the source began to transition to its flaring state, and the broadband fractional pulse amplitude increased above 7%, indicating an associated transition in the pulse waveform. The resulting pulse amplitude and phase spectra are shown in Figure 2.

Consistent with previous results (Patruno et al. 2009a; Sanna et al. 2017), we find that the pulse amplitude increases with energy up to ≈2 keV, and then flattens off (top panel). More interesting, however, is to note that the overall pulse energy spectrum shows significant evolution with time. As the source briefly rebrightens during its decay (second panel), the amplitude spectrum hardens, but the phase spectrum remains mostly the same. However, as the source intensity drops further, the phase residuals below 2 keV begin to pivot, while the fractional amplitudes increase again (third and bottom panels).

# 5. Long-term Evolution

## 5.1. Spin Frequency

Between 1998 and 2011, the spin frequency of SAX J1808 decreased at a constant rate (Patruno et al. 2012; for the 2015 outburst, see discussion). In Figure 3 we compare the three spin frequencies measured for the 2019 data with this established spin-down trend. The historic data for this comparison was adopted from Hartman et al. (2008, 2009) and Patruno et al.

(2012). The frequencies of the linear and flux-adjusted phase models both confirm a continuation of the long-term spin-down. The frequency derived from the quadratic phase model, however, lies well away from the trend. To quantitatively model the spin frequency evolution, we adopt the frequency from the linear phase model, as that model is the most similar to the historic data (see Section 6.1).

While there is a clear decreasing trend in the long-term evolution of the spin frequency, a linear function does not provide a statistically acceptable description of these data ($\chi_r^2$ of 5.1 for 5 dof). If we remove the long-term linear spin-down trend, and consider the residuals as a function of day of year, we observe an apparently sinusoidal modulation at the Earth's orbital period. This result suggests that our source coordinates are slightly offset from the optical position, such that the barycenter correction leaves a small periodic Doppler shift in the data (Manchester & Peters 1972; Manchester et al. 1974; Lyne & Graham-Smith 1990).

To account for the uncertainty in the source coordinates, we modeled the long-term evolution of the spin frequency as a linear trend plus a correction term for the residual Doppler shift,

$$ \Delta\nu(t) = \delta\nu_{98} + \dot{\nu}(t - T_{98}) + \mathcal{D}(t, \delta\lambda, \delta\beta), \quad (4) $$

where $\Delta\nu$ gives the measured spin frequency relative to $\nu_{98}$, the spin frequency measured for the epoch of the 1998 outburst, $T_{98}$. Additionally, $\delta\nu_{98}$ gives the correction to the spin frequency, $\dot{\nu}$ the spin frequency derivative at $T_{98}$, and $\mathcal{D}(t, \delta\lambda, \delta\beta)$ the Doppler correction term, which is defined as

$$ \mathcal{D}(t, \delta\lambda, \delta\beta) = \nu_p \frac{a_{\oplus}}{c} \omega [\delta\lambda \cos(\beta) \cos(\omega(t - T_{\Upsilon}) - \lambda) + \delta\beta \sin(\beta) \sin(\omega(t - T_{\Upsilon}) - \lambda)], \quad (5) $$

where $a_{\oplus}$ is the Earth's circular orbital radius, $c$ is the speed of light, $\omega = 2\pi/P_{\oplus}$ is the Earth's circular angular frequency, $P_{\oplus}$ is a sidereal year, and $T_{\Upsilon} = \text{MJD} 51079.143$ is the epoch of zero degrees ecliptic longitude for the Earth's orbit.<fnref n="12" /> Finally, $\lambda$ and $\beta$ give the respective source longitude and latitude in ecliptic coordinates, with $\delta\lambda$ and $\delta\beta$ their correction terms. We further note that corrections associated with the Earth's orbital eccentricity are negligible and can be safely ignored. Applying this model to the data, we obtain a good fit ($\chi_r^2$ of 1.7 for 3 dof) for a coordinate offset of $\delta\lambda = 0''33 \pm 0''10$ and $\delta\beta = -0''60 \pm 0''25$. In Figure 4 we show how this model fits the spin frequency residuals relative to the linear spin-down trend.

Our timing coordinates are offset from the optical position by a 2.3σ deviation. This offset is sufficiently large that one might wonder if the source shows measurable proper motion. To test for proper motion, we substituted the constant coordinate offsets in Equation (5) with their linear expansion in time, i.e.,

$$ \delta\lambda(t) = \delta\lambda_0 + \mu_{\lambda}(t - T_{\text{optical}}) \quad (6) $$

$$ \delta\beta(t) = \delta\beta_0 + \mu_{\beta}(t - T_{\text{optical}}), \quad (7) $$

where $T_{\text{optical}} = \text{MJD} 52703$ gives the 2001 June 13 date on which the optical data were collected (Hartman et al. 2008).

<footnote n="12">¹² This occurs when the Earth lies in the direction of the J2000 vernal equinox (\Upsilon) relative to the Sun, which is approximately the date of the northern autumnal equinox (September 23).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Energy resolved pulse fractional amplitude (left) and phase (right). Each panel shows a spectrum for the indicated time interval in black, with the other three intervals plotted in gray for reference.</figcaption></figure>

Assuming the optical position was exactly correct (i.e., fixing $\delta\lambda_0$ and $\delta\beta_0$ to zero), we obtain a very poor fit ($\chi_r^2$ of 7.1 for 3 dof). Letting all parameters vary yields a $\chi_r^2$ of 1.8 for 1 dof, with mean offset of $\delta\lambda_0 = 0''33 \pm 0''10$ and $\delta\beta_0 = -0''37 \pm 0''40$, and with ill-constrained derivatives of $\mu_\lambda = +0.013 \pm 0''10$ yr$^{-1}$ and $\mu_\beta = -0.04 \pm 0''03$ yr$^{-1}$. Hence, if the source shows a substantial proper motion, the X-ray data alone are insufficient to constrain it. Formally, the fixed position offset provides the best description of the data, so we transform that pulsar timing source position to equatorial coordinates, and list them in Table 2, along with the associated spin frequency and its derivative.

Finally, we note that similar results can be obtained by adopting the spin frequency obtained from the flux-adjusted phase model. In this case the long-term spin-down was slightly smaller, at $\dot{\nu} = (-8.5 \pm 0.5) \times 10^{-16}$ Hz s$^{-1}$, while the measured offset in ecliptic latitude reduced to $\delta\beta = -0''40 \pm 0''19$ and the offset in ecliptic longitude was $\delta\lambda = 0''28 \pm 0''07$.

## 5.2. Orbital Period

To study the long-term evolution of the orbital period, we follow the procedure of Hartman et al. (2008) and calculate the residual time of passage through the ascending node, $\Delta T_{asc}$, as

$$ \Delta T_{asc} = T_{asc,i} - (T_{ref} + NP_b), \qquad (8) $$

where we set the reference time to the epoch of the 2002 outburst, $T_{ref} = 52499.96$. Additionally, $T_{asc, i}$ is the time of ascending node measured for the $i$th outburst, $N$ is the number of orbital revolutions between the $i$th outburst and $T_{ref}$, and $P_b$ gives the orbital period reported by Hartman et al. (2009). These residual delays are shown in Figure 5. For context, we also show the parabolic model of Hartman et al. (2009), and the residuals with respect to this model.

The $\Delta T_{asc}$ measured for the 2019 outburst represents a significant lead with respect to the model of Hartman et al. (2009). From Figure 5 we see that this offset is a continuation of a recent trend, which suggests that the binary orbit has been contracting in the recent decade. Indeed, when we consider only the second decade of measurements (2008 through 2019),
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Spin frequency evolution of SAX J1808 as measured over two decades. Black bullets mark measurements made with RXTE (Hartman et al. 2008, 2009; Patruno et al. 2012). The colored squares give the NICER measurements for three phase models (see text). All frequencies are expressed relative to the 1998 epoch, with $\nu_{98} = 400.975210371$ Hz. The gray point marks the best measurement obtained for the 2015 outburst (Sanna et al. 2017); however, given its substantial systematic uncertainty (see the text), we did not include this point in our analysis. The solid line indicates the best-fit frequency evolution as determined in Section 5.1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Spin frequency measurements relative to the best-fit spin-down model (see also Figure 3) as a function of the Earth's ecliptic longitude. As the uncertainty of the 2015 outburst spin frequency span the vertical range, we replaced this point with a dashed line for clarity.</figcaption></figure>

<table><caption>**Table 2**
Long-term Timing Solution of SAX J1808</caption><thead><tr><td>Parameter</td><td>Value</td><td>Uncertainty</td></tr></thead><tbody><tr><td>R.A. (J2000)</td><td>18h08m27s647</td><td>0.007</td></tr><tr><td>decl. (J2000)</td><td>-36°58'43"90</td><td>0.025</td></tr><tr><td>Epoch (MJD)</td><td>50921.6</td><td>...</td></tr><tr><td>ν (Hz)</td><td>400.975210371</td><td>1.4 × 10<sup>-8</sup></td></tr><tr><td>ν (Hz s<sup>-1</sup>)</td><td>-1.01 × 10<sup>-15</sup></td><td>7 × 10<sup>-17</sup></td></tr></tbody></table>

**Note.** Uncertainties give the 1σ statistical errors. Both the spin frequency and its derivative are defined relative to the epoch of the 1998 outburst.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top: residual offset of the $T_{asc}$ measurements relative to predicted times based on the 2002 orbital period and reference time (Hartman et al. 2009). The curves show parabolic trends fit to the 1998–2008 (dashed), 2008–2019 (dashed–dotted), and 1998–2019 (solid) subsets of the data. Middle: residuals relative to the 1998–2008 parabolic model. Bottom: residuals relative to a parabolic model fit to the entire data set. A sinusoid with a 18.2 yr period and 7 s amplitude was added to guide the eye. In all panels the uncertainties are smaller than the point size.</figcaption></figure>

then we can fit the $\Delta T_{asc}$ evolution with a parabolic function ($\chi_r^2$ of 0.02 for 1 dof). This fit yields an orbital period derivative of $\dot{P}_b = (-5.18 \pm 0.02) \times 10^{-12} \text{ s s}^{-1}$.

Considering all data, we find that the mean long-term trend of the $\Delta T_{asc}$ yields a period derivative of $\dot{P}_b = (1.6 \pm 0.7) \times 10^{-12} \text{ s s}^{-1}$. The residuals around this mean trend show a large amplitude oscillation with a period of approximately 18.2 yr. This oscillation is clearly not sinusoidal, and could also not be adequately described as time delays due to a third body in an eccentric orbit. Instead, this quasi-periodic oscillation is apparently a stochastic process that is intrinsic to the dynamics of the binary orbit (see, e.g., Patruno et al. 2017; Sanna et al. 2017, for detailed modeling). What is not clear from the available data is whether this quasi-periodic process has a ≈7 samplitude around a steadily expanding orbit, or if we are seeing a ≈20 s amplitude modulation around a constant binary period. Additional monitoring of future outbursts is needed to differentiate between these two scenarios.

## 6. Discussion

We have presented a coherent timing analysis of the 2019 outburst of SAX J1808 as observed with NICER. We found that over the course of the outburst, the pulse phase changes as a function of the count rate. Considering the long-term evolution, we found there is now sufficient data for pulsar timing astrometry. In the following we consider some of the implications of these results.
## 6.1. Spin Frequency Evolution

We modeled the long-term evolution of the neutron star spin frequency by taking into account possible frequency shifts due to the position uncertainty. We found that this model gave a statistically acceptable description of the measurements, implying that we obtained an improved source position through pulsar timing astrometry. The timing position obtained through this analysis is consistent with the optical at a $2.3\sigma$ level. Additionally, we note that the timing position is consistent with the source coordinates in radio (Rupen et al. 2002) at a $1.4\sigma$ level. Even though the offset does not formally rise to a significant level, it is difficult to understand how the offset could be this large, given that the $0.15$ uncertainty in the optical position is dominated by the systematic uncertainty in the 2MASS catalog (Hartman et al. 2008). We investigated if the relatively large offset between the optical and pulsar timing positions could be explained through proper motion. However, we found that we could not constrain the magnitude of a potential proper motion from the pulsar timing data. Allowing for proper motion, however, did reduce the offset in coordinates. Because of the larger number of free parameters, it also increased the uncertainty on the pulsar timing position, which further reduced the tension between the various source position measurements. Hence, while a source position correction could explain the observed variations in long-term spin frequency evolution, this interpretation remains tentative. Future X-ray timing measurements or follow-up optical imaging are needed to test this interpretation.

An additional important caveat to this analysis is that the averaged spin frequency measured in each outburst is subject to systematic modeling uncertainty (as is clear from Table 1). For consistency, we only used spin frequency measurements that were obtained using a linear phase model (Hartman et al. 2008, 2009). While Patruno et al. (2009c) report slightly offset spin frequency measurements using a flux correction, there is an important distinction between the way we modeled the flux-adjustment and their work. Patruno et al. (2009c) determined the optimal timing solution by minimizing a linear fit to the phase-flux correlation, as a function of pulse frequency. In contrast, we incorporated the bias term directly into the phase model through a physically motivated power-law dependence. Hence, the reference frame of our flux-adjusted phase model is not guaranteed to be the same as that of Patruno et al. (2009c). Because a complete reanalysis of the archival RXTE data is beyond the scope of this work, we applied our spin-down analysis to the linear phase model. While this choice may have introduced some bias into our best-fit solution, the main implication of our analysis is that we now have a sufficiently long baseline of measurements for SAX J1808 that astrometric effects can no longer be treated as uncertainties, but instead should be accounted for in the timing analysis.

Finally, we note that we did not include a spin frequency measurement from the 2015 outburst. While a spin frequency was determined for this outburst using XMM-Newton and NuSTAR (Sanna et al. 2017), as well as with Chandra (Patruno et al. 2017), each of these measurements suffers from substantial statistical and systematic uncertainty. Specifically, the NuSTAR data suffered from clock drifts, which limited the absolute time resolution. The Chandra data, on the other hand, was collected late in the outburst, when the source had already decayed to low flux levels. These data have relatively poor statistical uncertainty on the frequency measurements. The

XMM-Newton data provided a smaller statistical uncertainty, but yielded very different spin frequencies when considering either the fundamental pulsation or its second harmonic (Sanna et al. 2017). While the measurement based on the second harmonic appears consistent with the long-term spin-down trend (Figure 3), the discrepancy between the two harmonics does suggest that the systematic uncertainty may be much larger than the statistical one. Additionally, all these observations only sample a very limited flux range, and are potentially biased by the flux-phase relation.

## 6.2. Hot-spot Drift

We found that, over the course of the 2019 outburst, the pulse phase showed significant evolution relative to a constant spin frequency model. We attempted to account for this evolution with two approaches: a quadratic phase model, and a flux-adjusted phase model.

The quadratic phase model accounts for the phase drift by adding a spin frequency derivative to the timing model. The physical interpretation of this model is that accretion torques applied during outburst are measurably affecting the spin frequency of the neutron star. While this model mostly accounts for the observed phase drift, there are significant issues with it. First, the residuals still show a (minor) systematic trend that is not consistent with timing noise. Second, the spin frequency derivative obtained through this model is anomalous. At $-3 \times 10^{-13} \text{ Hz s}^{-1}$, the derivative implies that the accretion torque is spinning-down the neutron star. Previous outbursts did not show such a spin-down. Although some early work on SAX J1808 report spin changes on the order of $10^{-13} \text{ Hz s}^{-1}$ (Morgan et al. 2003; Burderi et al. 2007), the sign and magnitude of the derivative change depending on how the data are selected, based on which Hartman et al. (2008) argued that the observed phase changes are actually rooted in the shape changes of the pulse waveform. Third, we expect that the accretion torque acts in the same manner for every outburst, meaning that each of the previous outbursts should have applied a similar spin-down torque to the neutron star. After accounting for the duration of each outburst, however, we find a long-term averaged accretion-powered spin-down of $-6 \times 10^{-15} \text{ Hz s}^{-1}$. This rate is significantly higher than the long-term averaged spin-down that is actually observed. Furthermore, the mean 2019 spin frequency obtained from this model is actually larger than the mean spin frequency measured during the 2011 outburst, further exacerbating this inconsistency. Given these contradictions, we conclude that the quadratic phase model is not physical.

In our second approach, we accounted for the phase shifts by adding a flux-dependent adjustment to our timing model. The physical interpretation for this phase model is that the location of the hot spot on the stellar surface is not fixed. Instead, the hot spot drifts as a function of mass accretion rate, causing apparent phase shifts that are not related to the stellar rotation (Romanova et al. 2003; Lamb et al. 2009; Patruno et al. 2009c). To describe the hot-spot drift, we adopted a power-law term in our phase model (Equation (3)). This flux-adjusted model performs as well as the quadratic model discussed above. The spin frequency obtained through this model is larger than that of the linear phase model (by 90 nHz, or $6\sigma$), but still broadly consistent with the long-term spin-down trend.

Our choice for a power-law dependence, and its $-1/5$ index, was guided by numerical simulations of accretion onto
magnetic stars (Kulkarni & Romanova 2013). These simulations indicate that the azimuthal position of the hot spot relative to the magnetic pole scales linearly with the magnetospheric radius as

$$ \phi_0 \propto -135^{\circ} \frac{r_m}{r_c}, \qquad (9) $$

where $r_m$ is the magnetospheric radius and $r_c$ the corotation radius. The latter is defined as the radius at which the Keplerian orbital frequency matches the spin frequency,

$$ r_c = \left( \frac{GM}{\Omega^2} \right)^{1/3}, \qquad (10) $$

with G the gravitational constant, M the neutron star mass, and $\Omega = 2\pi\nu_p$ the angular velocity of the neutron star. The former is defined as the inner truncation radius of the accretion disk, which follows from balancing the magnetospheric stresses against the material stresses in the accretion disk. The classical definition is (Pringle & Rees 1972; Lamb et al. 1973)

$$ r_m = \gamma_B^{2/7} \left( \frac{\mu^4}{GM\dot{M}^2} \right)^{1/7}, \qquad (11) $$

where $\mu$ is the magnetic dipole moment and $\gamma_B$ captures the uncertain physics determining the extent of the boundary layer where the disk couples to the magnetosphere (Psalis & Chakrabarty 1999). Numerical simulations, however, indicate a weaker dependence on mass accretion rate (Kulkarni & Romanova 2013),

$$ r_m = 1.06R \left( \frac{\mu^2}{GM\dot{M}R^7} \right)^{1/5}, \qquad (12) $$

with R the neutron star radius. The accretion rate can be estimated from the bolometric flux through the relation $4\pi d^2 F_{bol} = GM\dot{M}/R$. In Section 3 we further assumed that the bolometric flux, $F_{bol}$, is proportional to the observed count rate. To calibrate this relation, we fit each continuous good time interval using an absorbed multitemperature disk blackbody and a Comptonized power law (see also Papitto et al. 2009; Bult et al. 2019a; di Salvo et al. 2019). For each spectrum, we then estimated the bolometric flux by extrapolating the model between 0.01 keV and 100 keV. The resulting fluxes were indeed nearly proportional to the 0.5–10 keV count rate, $\Lambda$, with

$$ F_b = \kappa \Lambda, \qquad (13) $$

where we measured $\kappa = 4.14 \times 10^{-12} \text{ erg cm}^{-2} \text{ ct}^{-1}$. Rewriting these expressions, we find, for the azimuthal position of the hot spot,

$$ \phi_0 \propto -135^{\circ} \beta \Lambda^{-1/5}, \qquad (14) $$

where

$$ \beta = 2.92 \left[ \frac{M}{1.4 M_{\odot}} \right]^{-7/30} \left[ \frac{R}{10 \text{ km}} \right]^{1/10} \times \left[ \frac{\mu}{1.4 \times 10^{26} \text{ G cm}^3} \right]^{2/5} \left[ \frac{d}{3.5 \text{ kpc}} \right]^{-2/5}, \qquad (15) $$

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Neutron star mass–radius constraints from the bias scale of the flux-adjusted phase model. The dashed lines show isocontours for the magnetic dipole model parameter $\mu$ in units of $10^{26} \text{ G cm}^3$, assuming a distance of 3.5 kpc. The hatched region illustrates the distance uncertainty, for $\mu$ fixed at $1.0 \times 10^{26} \text{ G cm}^3$.</figcaption></figure>

Thus, from theory, we expect that the scale of our count-rate dependent bias should be $b = -1.09$, which is offset by about 20% from the $b = -0.87$ that we actually measured using our flux-adjusted phase model (Table 1). This offset can be accounted for by changing the assumed neutron star mass ($2M_{\odot}$), radius (8 km), and distance (3.8 kpc) (Galloway & Cumming 2006; Hartman et al. 2008; Goodwin et al. 2019a), but only if $\mu < 1.2 \times 10^{26} \text{ G cm}^3$ (see Figure 6).

The magnetic moment of SAX J1808 has been constrained to $(0.7 - 1.4) \times 10^{26} \text{ G cm}^3$ (Hartman et al. 2008; Patrino et al. 2012). The upper limit follows from relating the long-term spin-down trend to magnetic dipole radiation. The lower limit is obtained by requiring that, at the highest luminosities, the magnetic field be strong enough to truncate the disk at the neutron star radius (Psalis & Chakrabarty 1999), although we note this value changes depending on specific assumptions for the source distance and the neutron star parameters (see, e.g., Mukherjee et al. 2015, for an overview). In this context, the upper limit on $\mu$ that we obtained through the bias scale may indicate that the magnetic dipole radiation is not the only source of torque applied to the neutron star. Instead, an additional torque may be present in the form of the propeller mechanism (Illarionov & Sunyaev 1975; D’Angelo & Spruit 2010), or gravitational wave emission due to a neutron star mass quadrupole (Wagoner 1984; Bildsten 1998; Anderson & Kokkotas 2001; Melatos & Payne 2005; Mahmoodifar & Strohmayer 2013).

Alternatively, the tension between the magnetic moment from the spin-down trend and that from the bias scale could be attributed to modeling uncertainties related to the disk–magnetosphere boundary region. By determining the point at which the magnetic field is strong enough to force the disk into corotation, D’Angelo & Spruit (2010) derive the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Neutron star mass–radius constraints from the bias scale of the flux-adjusted phase model. The dashed lines show isocontours for the magnetic diffusivity parameter η, assuming a distance of 3.5 kpc. The hatched region illustrates the distance uncertainty, for η fixed at 0.3. The magnetic moment has been fixed to μ = 1.4 × 10²⁶ G cm³.</figcaption></figure>

magnetospheric radius as

$$r_m = \left( \frac{\eta \mu^2}{4 \Omega \dot{M}} \right)^{1/5}, \qquad (16)$$

where η parameterizes the poorly understood magnetic diffusivity of the accretion disk. We point out that this relation almost exactly matches Equation (12) when η = 1 and the neutron star mass and radius are fixed to the values used in the simulations of Kulkarni & Romanova (2013). The different dependence on R is due to the assumption of what the orbital velocity at the inner edge of the disk should be. Using Equation (16) we can derive an alternative expression for β as

$$\beta = 2.92 \times \eta^{1/5} \left[ \frac{M}{1.4 M_{\odot}} \right]^{-2/5} \left[ \frac{R}{10 \text{ km}} \right]^{-1/5} \times \left[ \frac{\mu}{1.4 \times 10^{26} \text{ G cm}^3} \right]^{2/5} \left[ \frac{d}{3.5 \text{ kpc}} \right]^{-2/5} . \quad (17)$$

Assuming magnetic dipole radiation is the dominant source of spin-down torque (i.e., setting μ = 1.4 × 10²⁶ G cm³), we can constrain η to ≈0.3–0.4 (see Figure 7).

While these constraints are subject to the various uncertainties in the numerical modeling, the fact that we can successfully describe the systematic timing noise using a physical model provides a strong argument in favor of this interpretation. Additionally, that we can obtain reasonable values for the properties of the neutron star bodes well for future efforts in detailed waveform modeling of accreting millisecond pulsars. In this analysis we only consider the phase residuals, but additional information is encoded in the precise shape of the waveform (Poutanen & Gierliński 2003; Ibragimov & Poutanen 2009; Morsink & Leahy 2011), and in its energy dependence. Hence, comprehensive time-resolved waveform

modeling of SAX J1808 may not only reveal the neutron star mass and radius, but also constrain the detailed physics of the disk–magnetosphere boundary region.

This work was supported by NASA through the NICER mission and the Astrophysics Explorers Program, and made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC). P.B. was supported by an NPP fellowship at NASA Goddard Space Flight Center.

Facilities: ADS, HEASARC, NICER.

Software: heasoft (v6.26), nicerdas (v6a), tempo2 (Hobbs et al. 2006).

## ORCID iDs

Deepto Chakrabarty https://orcid.org/0000-0001-8804-8946

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

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

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

Tod E. Strohmayer https://orcid.org/0000-0001-7681-5845

## References

Andersson, N., & Kokkotas, K. D. 2001, *IJMPD*, 10, 381

Arzoumanian, Z., Fruchter, A. S., & Taylor, J. H. 1994, *ApJL*, 426, L85

Bhattacharyya, S., & Strohmayer, T. E. 2007, *ApJ*, 656, 414

Bildsten, L. 1998, *ApJL*, 501, L89

Bult, P., Jaisawal, G. K., Güver, T., et al. 2019a, *ApJL*, 885, L1

Bult, P., & van der Klis, M. 2015a, *ApJL*, 798, L29

Bult, P., & van der Klis, M. 2015b, *ApJ*, 806, 90

Bult, P. M., Gendreau, K. C., Arzoumanian, Z., et al. 2019b, *ATel*, 13077, 1

Bult, P. M., Gendreau, K. C., Arzoumanian, Z., et al. 2019c, *ATel*, 13001, 1

Burderi, L., di Salvo, T., Lavagetto, G., et al. 2007, *ApJ*, 657, 961

Burderi, L., Riggio, A., di Salvo, T., et al. 2009, *A&A*, 496, L17

Cackett, E. M., Altamirano, D., Patruno, A., et al. 2009, *ApJL*, 694, L21

Campana, S., & di Salvo, T. 2018, in Astrophysics and Space Science Library, Vol. 457, The Physics and Astrophysics of Neutron Stars, ed. L. Rezzolla et al. (Cham: Springer), 149

Chakrabarty, D., & Morgan, E. H. 1998, *Natur*, 394, 346

Chakrabarty, D., Morgan, E. H., Muno, M. P., et al. 2003, *Natur*, 424, 42

D'Angelo, C. R., & Spruit, H. C. 2010, *MNRAS*, 406, 1208

di Salvo, T., Burderi, L., Riggio, A., Papitto, A., & Menna, M. T. 2008, *MNRAS*, 389, 1851

di Salvo, T., Sanna, A., Burderi, L., et al. 2019, *MNRAS*, 483, 767

Doroshenko, O., Löhmer, O., Kramer, M., et al. 2001, *A&A*, 379, 579

Galloway, D. K., & Cumming, A. 2006, *ApJ*, 652, 559

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

Gierliński, M., Done, C., & Barret, D. 2002, *MNRAS*, 331, 141

Goodwin, A. J., Galloway, D. K., Heger, A., Cumming, A., & Johnston, Z. 2019a, *MNRAS*, 490, 2228

Goodwin, A. J., Russell, D. M., Galloway, D. K., et al. 2019b, *ATel*, 12993, 1

Hartman, J. M., Patruno, A., Chakrabarty, D., et al. 2008, *ApJ*, 675, 1468

Hartman, J. M., Patruno, A., Chakrabarty, D., et al. 2009, *ApJ*, 702, 1673

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

Ibragimov, A., & Poutanen, J. 2009, *MNRAS*, 400, 492

Illarionov, A. F., & Sunyaev, R. A. 1975, *A&A*, 39, 185

in't Zand, J. J. M., Galloway, D. K., Marshall, H. L., et al. 2013, *A&A*, 553, A83

in't Zand, J. J. M., Heise, J., Muller, J. M., et al. 1998, *A&A*, 331, L25

in't Zand, J. J. M., van Kerkwijk, M. H., Pooley, D., et al. 2001, *ApJL*, 563, L41

Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, *ApJS*, 163, 401

Kulkarni, A. K., & Romanova, M. M. 2013, *MNRAS*, 433, 3048

Lamb, F. K., Boutloukos, S., van Wassenhove, S., et al. 2009, *ApJL*, 705, L36

Lamb, F. K., Pethick, C. J., & Pines, D. 1973, *ApJ*, 184, 271

Lyne, A. G., & Graham-Smith, F. 1990, Pulsar Astronomy (Cambridge: Cambridge Univ. Press)

Mahmoodifar, S., & Strohmayer, T. 2013, *ApJ*, 773, 140

Manchester, R. N., & Peters, W. L. 1972, *ApJ*, 173, 221

Manchester, R. N., Taylor, J. H., & Van, Y. Y. 1974, *ApJL*, 189, L119
Melatos, A., & Payne, D. J. B. 2005, ApJ, 623, 1044

Morgan, E. H., Chakrabarty, D., Wijnands, R., van der Klis, M., & Markwardt, C. 2003, AAS/High Energy Astrophysics Division Meeting, 7, 17.29

Morsink, S. M., & Leahy, D. A. 2011, ApJ, 726, 56

Mukherjee, D., Bult, P., van der Klis, M., & Bhattacharya, D. 2015, MNRAS, 452, 3994

Papitto, A., di Salvo, T., D'Aì, A., et al. 2009, A&A, 493, L39

Parikh, A. S., & Wijnands, R. 2019, ATel, 13000, 1

Patruno, A., Bult, P., Gopakumar, A., et al. 2012, ApJL, 746, L27

Patruno, A., Jaodand, A., Kuiper, L., et al. 2017, ApJ, 841, 98

Patruno, A., Maitra, D., Curran, P. A., et al. 2016, ApJ, 817, 100

Patruno, A., Rea, N., Altamirano, D., et al. 2009a, MNRAS, 396, L51

Patruno, A., Watts, A., Klein Wolt, M., Wijnands, R., & van der Klis, M. 2009b, ApJ, 707, 1296

Patruno, A., & Watts, A. L. 2012, arXiv:1206.2727

Patruno, A., Wijnands, R., & van der Klis, M. 2009c, ApJL, 698, L60

Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301

Pringle, J. E., & Rees, M. J. 1972, A&A, 21, 1

Psaltis, D., & Chakrabarty, D. 1999, ApJ, 521, 332

Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., Wick, J. V., & Lovelace, R. V. E. 2003, ApJ, 595, 1009

Rupen, M. P., Dhawan, V., Mioduszewski, A. J., Stappers, B. W., & Gaensler, B. M. 2002, IAUC, 7997, 2

Russell, D. M., Bramich, D. M., Lewis, F., et al. 2019a, AN, 340, 278

Russell, D. M., Goodwin, A. J., Galloway, D. K., et al. 2019b, ATel, 12964, 1

Sanna, A., di Salvo, T., Burderi, L., et al. 2017, MNRAS, 471, 463

Standish, E. M. 1998, JPL Planetary and Lunar Ephemerides, DE405/LE405, JPL Interoffice Memo 312.F-98-048 (Pasadena, CA: NASA Jet Propulsion Laboratory)

Wagoner, R. V. 1984, ApJ, 278, 345

Wijnands, R., Méndez, M., Markwardt, C., et al. 2001, ApJ, 560, 892

Wijnands, R., & van der Klis, M. 1998, Natur, 394, 344