{'fraction_non_alphanumeric': 0.03582897240723121, 'fraction_numerical': 0.03279614652711703, 'mean_word_length': 4.555500495540139, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Timing the pulsations of the accreting millisecond pulsar SAX J1808.4-3658 during its 2019 outburst

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, enabling pulsar timing based astrometry of this accreting millisecond pulsar.


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, semi-regular outburst recurrence times of 2 − 4 years. 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,b), its spectral characteristics (Gierliński * NASA Postdoctoral Fellow 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(Hartman et al. , 2009Patruno et al. , 2017Sanna et al. 2017), while between outbursts the source spins down at a steady rate of ≈ 1.5 × 10 −15 Hz s −1 . Additionally, these analyses revealed the 2.01 hr binary orbit (Chakrabarty & Morgan 1998), and estab-lished 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 −12 s s −1 (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 spinorbit angular momentum exchange in the companion star (Hartman et al. 2008;), 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 and2015, showed that the rate of orbital expansion is itself evolving rapidly. These measurements have been analysed 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. (2019a) 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 . 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.

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. 2009bPatruno et al. , 2016Bult 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 inter-vals where the rate of saturating particle events (overshoots) is greater than 1 ct s −1 detector −1 , or greater than 1.52 × cor sax −0.633 , where cor sax 1 gives the cut-off rigidity of the Earth's magnetic field, in units of GeV c −1 . 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 cut-off threshold for the overshoots to 1.5 ct s −1 detector −1 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 (TDB). 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.

COHERENT TIMING
We performed a coherent timing analysis of the 401 Hz pulsations of SAX J1808. For each continuous goodtime 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 ν k = kν p , where ν p is the pulsar spin frequency. Additionally, for each sinusoid we determined the phase, ϕ k , and amplitude, A k , from which we derive a fiducial pulse arrival time and a fractional sinusoidal amplitude, r, as 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, 1st and 2nd 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, and 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 semi-major axis, a x sin i, the time of passage through the ascending node, T asc , and the spin frequency, ν 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 χ 2 (χ 2 r ) 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, 3rd panel). This trend suggests that the spin frequency is changing over the course of the outburst. Indeed, including a spin frequency derivative in our timing model improves the fit statistic to a χ 2 r 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, 4th 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 where ϕ 0 gives a reference phase, ν p t gives the linear term due to the constant spin frequency, ϕ orb (t) gives the phase correction associated with the binary orbit, and ϕ bias (F bol ) gives a phase-bias correction term as a function of bolometric flux, F 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 aṡ M −1/5 (D'Angelo & Spruit 2010; Kulkarni & Romanova 2013). SinceṀ should be proportional to the bolometric flux, we define a simple power-law model with Γ 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 χ 2 r of 2.08 for 258 dof. While this fit statistic is still not formally acceptable, it performs better than the quadratic phase model. Furthermore, we note that phase residuals for this model (Figure 1, 5th panel) are no longer structured, such that the poor statistic can be attributed to residual timing noise (Patruno & Watts 2012).
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.

PULSE ENERGY DEPENDENCE
To investigate the photon-energy dependence of the pulsations, we applied a sliding window method in energy space. That is, we selected all data in a 0.25 keV wide window and measured the resulting pulse profile properties. We then repeated this analysis by moving the window from 0.5 keV to 10 keV using steps of 0.1 keV. This method was applied to four representative time intervals. The first interval covers August 10 through 17, which is when the outburst peaked and the broad-band fractional pulse amplitude remained constant at 5%. The second interval covers August 18 through 20, which is where the source flux showed a re-brightening with respect to its long-term decay, and, likewise, the broad-band 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 broad-band fractional pulse amplitude increased above 7%, indicating an associated tran- Flux-adjusted phase model Note-Uncertainties give the 1σ statistical errors and the upper limit is quoted at the 95% c.l. The parameter gives the orbital eccentricity.
sition 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 re-brightens during its decay (second panel), the amplitude spectrum hardens, but the phase remains mostly the same. However, as the source intensity drops further, the phase residuals below 2 keV begin to pivot, as the fractional amplitudes below 2 keV show large variations (third and bottom panels).

LONG TERM EVOLUTION 5.1. Spin frequency
Between 1998 and 2011, the spin frequency of SAX J1808 decreased at a constant rate ; 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. (2008Hartman et al. ( , 2009) and . 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 (χ 2 r of 5.1 for 5 dof). If we remove the longterm 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, where ∆ν gives the measured spin frequency relative to ν 98 , the spin frequency measured for the epoch of the 1998 outburst, T 98 . Additionally, δν 98 gives the correction to the spin frequency,ν the spin frequency derivative at T 98 , and D(t, δλ, δβ) the Doppler correction term, which is defined as where a ⊕ is the Earth's circular orbital radius, c is the speed of light, ω = 2π/P ⊕ is the Earth's circular angular frequency, P ⊕ is a sidereal year, and T Υ = MJD 51079.143 is the epoch of zero degrees ecliptic longitude for the Earth's orbit 2 . Finally, λ and β give the respective source longitude and latitude in ecliptic coordinates, with δλ and δβ 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 (χ 2 r of 1.7 for 3 dof) for a coordinate offset of δλ = 0.33 ±0.10 and δβ = −0.60 ±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 2 This occurs when the Earth lies in the direction of the J2000 vernal equinox (Υ) relative to the Sun, which is approximately the date of the northern autumnal equinox (September 23). Year of outburst where T optical = MJD 52703 gives the 2001 June 13 date on which the optical data were collected (Hartman et al. 2008). Assuming the optical position was exactly correct (i.e., fixing δλ 0 and δβ 0 to zero), we obtain a very poor fit (χ 2 r of 7.1 for 3 dof). Letting all parameters vary yields a χ 2 r of 1.8 for 1 dof, with mean offset of δλ 0 = 0.33 ± 0.10 and δβ 0 = −0.37 ± 0.40 , and with ill-constrained derivatives of µ λ = +0.013 ± 0.10 yr −1 and µ β = −0.04 ± 0.03 yr −1 . Hence, the long-term spin frequency measurements are best described using a fixed position offset. We transform the 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 fluxadjusted phase model. In this case the long term spin down was slightly smaller, atν = (−8.5 ± 0.5) × 10 −16 Hz s −1 , while the measured offset in ecliptic latitude reduced to δβ = −0.40 ± 0.19 and the offset in ecliptic longitude was δλ = 0.28 ± 0.07 .   Note-Uncertainties give the 1σ statistical errors. Both the spin frequency and its derivative at are defined relative to the epoch of the 1998 outburst.

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, ∆T asc , as 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 ith outburst, N is the number of orbital revolutions between the ith 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 ∆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), then we can fit the ∆T asc evolution with a parabolic function (χ 2 r of 0.02 for 1 dof). This fit yields an orbital period derivative ofṖ b = (−5.18 ± 0.02) × 10 −12 s s −1 .
Considering all data, we find that the mean long-term trend of the ∆T asc yields a period derivative ofṖ b = (1.6±0.7)×10 −12 s s −1 . The residuals around this mean trend show a large amplitude oscillation with a period of approximately 18.2 years. 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 s amplitude 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.

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.

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σ level. Given that the 0.15 uncertainty in the optical position is set 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 no significant proper motion could be measured from the pulsar timing data.
An important caveat to this analysis is that the averaged spin frequency measured in each outburst is subject to systematic modelling 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(Hartman et al. , 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 fluxadjustment 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 spindown 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 XMM-Newton and 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. Additionally, all these observations only sample a very limited flux range, and are potentially biased by the flux-phase relation.

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 × 10 −13 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 Hz s −1 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, so by rescaling the apparent spin-down measured in 2019 based on the outburst recurrence cycle, we obtain a long-time averaged accretion-powered spin down of −6 × 10 −15 Hz s −1 . This rate is significantly higher than the long-term 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;  2009c). To describe the hot-spot drift, we adopted a power-law term in our phase model (Equation 3). This flux-adjusted model performs significantly better than 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σ), 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 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, with G the gravitational constant, M the neutron star mass, and Ω = 2πν p the angular spin frequency 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) where µ is the magnetic dipole moment and γ B captures the uncertain physics determining the extent of the boundary layer where the disk couples to the magnetosphere (Psaltis & Chakrabarty 1999). Numerical simulations, however, indicate a weaker dependence on mass accretion rate (Kulkarni & Romanova 2013), with R the neutron star radius. The accretion rate can be estimated from the bolometric flux through the relation 4πd 2 F bol = GMṀ /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 multi-temperature disk blackbody and a Compontized power-law (see also, Papitto et al. 2009;Di Salvo et al. 2019;Bult et al. 2019a). 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 countrate, Λ, with where we measured κ = 4.14 × 10 −12 erg cm −2 ct −1 .
Rewriting these expressions, we find, for the azimuthal position of the hot-spot, Thus, from theory, we expect that the scale of our countrate 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 by varying the neutron star parameters within their allowed ranges (Galloway & Cumming 2006;Hartman et al. 2008;Goodwin et al. 2019a), but only if µ < 1.2 × 10 26 G cm 3 (see Figure 6). The magnetic moment of SAX J1808 has been constrained to (0.7−1.4)×10 26 G cm 3 (Hartman et al. 2008;. 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 (Psaltis & 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 µ 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;Andersson & 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 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 co-rotation, D'Angelo & Spruit (2010) derive the magnetospheric radius as where η parameterizes the poorly understood magnetic diffusivity of the accretion disk. We point out that this relation almost exactly matches Equation 11 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 15 Assuming magnetic dipole radiation is the dominant source of spin-down torque (i.e., setting µ = 1.4 × 10 26 G cm 3 ), we can constrain η to ≈ 0.3−0.4 (see Figure  7). While these constraints are subject to the various uncertainties in the numerical modelling, 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 modelling 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 modelling 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.