Timing Analysis of the 2022 Outburst of SAX J1808.4-3658: Hints of Orbital Decay

We present a pulse timing analysis of NICER observations of the accreting millisecond X-ray pulsar SAX J1808.4 − 3658 during the outburst that started on 2022 August 19. Similar to previous outbursts, after decaying from a peak luminosity of ; 1 × 10 36 erg s − 1 in about a week, the pulsar entered a ∼ 1 month long re ﬂ aring stage. Comparison of the average pulsar spin frequency during the outburst with those previously measured con ﬁ rmed the long-term spin derivative of


Introduction
The transient low-mass X-ray binary SAX J1808.4−3658 (hereafter SAX J1808) was discovered in 1996 with the X-ray satellite BeppoSAX during an X-ray outburst (Zand et al. 1998). Three type I X-ray bursts were detected (in 't Zand et al. 2001), permitting the identification of the accretor as a neutron star (NS). The distance was later estimated by Galloway & Cumming (2006) to be ∼3.5 kpc. The detection of 401 Hz X-ray pulsations with RXTE during the following outburst in 1998 marked the discovery of the first accreting millisecond X-ray pulsar (AMXP; Wijnands & van der Klis 1998). Timing analysis revealed that the NS is in an orbit with a ≈0.05 M e brown dwarf companion (Bildsten & Chakrabarty 2001) with a 2.01 hr orbital period (Chakrabarty & Morgan 1998). Since its discovery, the source has undergone ten ∼1 month long outbursts with ∼2-3 yr recurrence. This makes it the AMXP that has shown the largest number of outbursts of sufficient duration for in-depth investigation of its long-term timing properties (Marino et al. 2019;Di Salvo & Sanna 2022). It is thus the most thoroughly studied AMXP. The X-ray luminosity typically reaches L X ∼ few × 10 36 erg s −1 (Gilfanov et al. 1998) at the peak of the outbursts, and decreases down to L X ∼ few × 10 31 erg s −1 during quiescence (Stella et al. 2000;Campana et al. 2004). Coherent 401 Hz X-ray pulsations are observed only during the outbursts and interpreted in terms of magnetic channeling of the in-flowing matter onto the NS magnetic poles. During the 2019 outburst, Ambrosino et al. (2021) discovered coherent millisecond optical and UV pulsations. The bright pulsed luminosity (L optical ≈ 3 × 10 31 erg s −1 , and L UV ≈ 2 × 10 32 erg s −1 ) seen at those wavelengths challenged the expectations of the standard accretion models.
On 2022 August 19, the MAXI/GSC nova alert system (Imai et al. 2022) detected the occurrence of a new outburst of SAX J1808, later confirmed by rapid targeted follow-up NICER observations (Sanna et al. 2022a). Here, we report on the high-cadence monitoring campaign performed with the NICER guest observer program ID 5574 (PI: A. Papitto). We focus on the pulse phase timing analysis carried out on the X-ray pulsations detected throughout the outburst to measure the pulsar spin frequency and the binary orbital ephemeris. These values are compared with those observed during previous outbursts (Di Salvo et al. 2008;Hartman et al. 2008;Burderi et al. 2009;Hartman et al. 2009;Patruno et al. 2016;Sanna et al. 2017;Bult et al. 2020) to derive the longterm spin and orbital evolution of the pulsar and discuss the implications for the models of AMXPs.
The NICER monitoring started when SAX J1808 had almost attained its peak count rate of ∼300 c s −1 (top panel of Figure 1).
To estimate the peak luminosity, we extracted the spectrum collected in the observation on 2022 August 24 (ObsID 5574010102), and modeled it within the XSPEC spectral fitting package (Arnaud 1996). We accounted for absorption effects using the tbabs model with wilm abundances (Wilms et al. 2000) and vern cross sections (Verner et al. 1996). We described the continuum emission using a combination of a disk blackbody (diskbb) and a Comptonization component (nthComp), adding three Gaussian emission lines. The electron temperature was held fixed to 30 keV in the fit. We obtained a satisfactory fit (χ 2 /dof = 865.70/850). We then calculated the 0.6-10 keV X-ray unabsorbed flux using the convolution model cflux. The corresponding peak luminosity is ≈1 × 10 36 erg s −1 (assuming a distance of 3.5 kpc; Galloway & Cumming 2006).
After ∼5 days from the first detection, the decay phase begun, until the source entered its typical reflaring stage (Cui et al. 1998;Wijnands et al. 2001;Hartman et al. 2008;Patruno & Watts 2021), which was observed with NICER for more than a month (Illiano et al. 2022). The light curve of the 2022 outburst slightly differed from the typical profile shown by SAX J1808. The usually shortlived peak exhibited the longest duration observed so far, and the slow decay/rapid drop lasted much less (∼10-15 days) than usual. No type I X-ray bursts were detected during NICER observations, unlike most other outbursts (see, e.g., in't Zand et al. 2001;Chakrabarty et al. 2003;Galloway & Cumming 2006;Bult & Klis 2015;Bult et al. 2020).

Coherent Timing
In order to correct the photon arrival times for the pulsar orbital motion in the binary system, we first performed a preliminary search on the epoch of passage at the ascending node, T asc , exploiting the variance of the epoch-folding search as a statistical estimator. We fixed the orbital period and the projected semimajor axis equal to the values found in the timing solution of the 2019 outburst (Bult et al. 2020), and we used the best T asc found as a starting point for the pulse phase timing. After correcting the photon arrival times with this preliminary orbital solution, we divided our data set into 1000 s Top panel: the 0.5-10 keV light curve using 200 s bins. Second panel: pulse fractional amplitude for the first harmonic (black) and the second harmonic (red). Third and fourth panels: phase residuals relative to the linear phase model for the first harmonic (black) and the second harmonic (red), and flux-adjusted phase models, respectively (also see Table 1). The phase residuals relative to the quadratic model are not plotted as they are similar to those of the linear model (see text). long segments and folded them around our best estimate of the spin frequency ν F using 16 phase bins. We modeled our pulse profiles with a constant plus two harmonic components, retaining only data in which the signal was detected with an amplitude significant at more than a 3σ confidence level. Throughout the outburst, the amplitude of the fundamental (black dots in the second panel of Figure 1) is higher than the second harmonic (red dots in the same panel), increasing by approximately ∼1-2 percentage points when the rapid drop phase of the outburst took place and slightly again at the beginning of the flaring tail.
We modeled the time evolution of the phase of the fundamental, using (see, e.g. Burderi et al. 2007;Papitto et al. 2007;Sanna et al. 2022b): Here, ν is the pulsar spin frequency, T 0 is the chosen reference epoch, f 0 is the pulse phase at T 0 , Δν = ν(T 0 ) − ν F , while R orb (t) is the residual Doppler modulation due to a difference between the adopted orbital parameters and the actual ones (see, e.g., Deeter et al. 1981). Table 1 shows the best-fitting orbital and spin parameters we obtained. To take into account the large value of the reduced χ 2 obtained from the fit, we rescaled the uncertainties of the fit parameters by the square root of that value (see, e.g., Finger et al. 1999). We estimated the systematic uncertainty on the spin frequency due to the positional uncertainties of the source using the expression ( ) y P 1 sin 2 where y = r E /c is the semimajor axis of the Earth orbit in light seconds, σ γ is the positional error circle, β is the source latitude in ecliptic coordinates, and P ⊕ is the Earth orbital period (see, e.g., Lyne & Graham-Smith 1990;Burderi et al. 2007;Sanna et al. 2017). Adopting the positional uncertainties reported by Bult et al. (2020), we estimate 5 10 Hz 8 pos sń - . We added in quadrature this systematic uncertainty to the statistical error of the spin frequency reported in Table 1.
We base the discussion of the phase evolution during the 2022 outburst on the properties of the fundamental frequency.
In fact, below 3 keV, the second harmonic was often too weak to be detected by NICER Bult et al. 2020). We modeled the phase delays using either a constant frequency model (i.e., setting  0 n = in Equation (1); see Table 1) or a constant spin frequency derivative. The quadratic fit returns a value of the average frequency derivative of  n = 2.4(4.0) × 10 −15 Hz s −1 (χ 2 /dof = 698.2/284), which is compatible with zero. The probability of a chance improvement of the χ 2 compared to the constant frequency model obtained with an F-test is ∼0.5, indicating that the addition of such a component does not produce a significant improvement in the data description.
A strong variability of the phase and shape of the pulse profiles characterized all SAX J1808 outbursts observed so far (see the reviews by Patruno & Watts 2021;Di Salvo & Sanna 2022). This strongly limited the ability to measure the NS spin evolution during individual outbursts from pulse phase timing. Pulse phases measured from the second harmonic generally showed a more regular behavior compared to the fundamental. Burderi et al. (2006) exploited this property to infer a spin-up rate of ( )  4.4 8 10 13 n =´-Hz s −1 during the 2002 outburst. Such a value is only slightly larger than that expected considering the material torque exerted by accretion through a Keplerian disk in-flow truncated a few tens of kilometers from the NS (;2 × 10 −13 Hz s −1 for a 1.4 M e NS accreting at a rate of 10 −9 M e yr −1 from a disk truncated at 20 km from the NS; see, e.g., Di Salvo et al. 2019 and references therein). Hartman et al. (2008Hartman et al. ( , 2009) attributed instead much of the observed phase variability to a red noise process affecting the pulse phases on timescales similar to the outburst duration; this led to tighter upper limits on the spin frequency derivative (| |  2.5 10 14 n <´-Hz s −1 ). Patruno et al. (2009) characterized such a noise process in terms of a correlation between the pulse phase and the X-ray flux. Azimuthal drifts of the hot spot location on the NS surface related to a movement of the inner disk truncation radius at different mass accretion rates could explain such a phase-flux correlation. However, a broadly different correlation characterized each of the outbursts of SAX J1808. In this context, Bult et al. (2020) found the best description of the evolution of the pulse phases measured by NICER in 2019 by using a phase-flux correlation term related to hot spots drifts, (3), and Γ = −0.2 fixed. These values were broadly consistent with those expected according to numerical simulations of accretion onto a fast-rotating NS. The fixed power-law index arises from the linear scaling of the azimuthal position of the hot spot with the magnetospheric radius, which was recently predicted to depend on the mass accretion rate as  M 1 5 - (Kulkarni & Romanova 2013). Because of the large phase residuals with respect to the linear model, following Bult et al. (2020) we also attempted to replace in Equation (1) the spin frequency derivative term with a component including a dependence of the pulse phase on the flux, here considered to be traced by the 0.5-10 keV count rate . The resulting χ 2 shows a significant improvement with respect to the linear phase model (F-test probability of ∼8.5 × 10 −28 ; see also panel 4 in Figure 1).
If we restrict our analysis to the first ∼8 days of the outburst, i.e., until the source faded to roughly a fifth of the peak flux, then the addition of either a spin frequency derivative or of a phase-flux correlation component did not improve the phase description compared to a constant spin frequency model (F-test probability of the quadratic model with respect to the linear one of ∼0.7; F-test probability of the flux-adjusted model with respect to the linear one of ∼0.8). The phase behavior is compatible with a constant spin frequency, with a 90% C.L. upper limit on the spin frequency derivative of 1.9 × 10 −13 Hz s −1 (same order of magnitude as the expected one for the accretion-driven spin-up, as discussed above). This points to an anticorrelation between the phase delays and the source flux, observed in Figure 1, holding only for count rates lower than ∼100 c s −1 , i.e., in the reflaring phase. Even though we lacked a coverage of the rising part of the outburst, i.e., where most of the flux dependence of the phases was present in 2019 data (see Figure 1 in Bult et al. 2020), we found an even more pronounced phase-flux anticorrelation than in the previous outburst.
Since the Γ index we obtained is not consistent with the hot spots drifts predicted by numerical simulations of accreting pulsars (Kulkarni & Romanova 2013), the phase shifts are not driven by the changing size of the magnetosphere, but are instead inversely proportional to the mass accretion rate (similar to the case of the AMXP MAXI J1816−195; Bult et al. 2022). On the other hand, no such variation was seen when the flux varied by a three-times larger factor during the peak and the decay phase. The steep index of the phase-flux correlation we found (δf ∼ 1/F X ) naturally explains why introducing this term determines a significant improvement of the quality of the residuals of the fit performed on the whole data set, even though phase fluctuations are essentially observed only at low count rates.

Long-term Spin Frequency Evolution
The ten SAX J1808 outbursts observed so far, the most numerous for any AMXP, enable a detailed study of the longterm spin frequency evolution through a comparison of the measurements obtained in each of the outbursts. Previous works (see, e.g., Patruno et al. 2012;Sanna et al. 2017;Bult et al. 2020) found that the spin frequency decreased at an average rate of   10 SD 15 n --Hz s −1 , compatible with the energy losses expected from a ≈10 26 G cm 3 rotating magnetic dipole. Bult et al. (2020) also found that the spin frequencies measured by correcting the pulse arrival times with the position measured by Hartman et al. (2008) showed a yearly modulation due to an offset of δλ = (0.33 ± 0.10)″ and δβ = (−0.60 ± 0.25)″ in the assumed Galactic longitude and latitude of the pulsar, respectively. In order to compare the frequency observed in the 2022 outburst with the past values, in this part of the analysis we corrected the photon arrival times to the SSB adopting the optical coordinates from Hartman et al. (2008), in analogy with previous works. Using a linear phase model, we obtained ν = 400.9752095863(45) Hz, higher than ∼8 × 10 −7 Hz compared to the values obtained with the coordinates from Bult et al. (2020). We then modeled the long-term frequency evolution (see Figure 2) with a function including a constant spin-down and a position correction term:

Orbital Period Evolution
To investigate the orbital evolution, we computed the difference   (Bult et al. 2020). The red dot is from this work, having corrected the data with the source coordinates from Hartman et al. (2008) and fitted the phase delays with a linear model. All frequencies are expressed relative to the 1998 spin frequency, ν 98 = 400.975210371 Hz (Hartman et al. 2008). The dotted line indicates the best-fitting function including the Doppler modulation due to source coordinates error, and the dashed line is the corresponding linear function. Bottom panel: residuals relative to the best-fitting function. We did not include the 2015 spin frequency estimate because its uncertainty is about a factor of 20 larger than the others and compatible with the amplitude of Doppler modulation. a constant period derivative (dotted line in Figure 3), leaves evident residuals with a sinusoidal shape (χ 2 /dof =15579.0/6, see the middle panel of Figure 3). We then added a sinusoidal term to the relation used to fit the orbital phases: 1 --, A = 11.30(33) s, and P = 7.57(21) × 10 3 days. We evaluated the uncertainties by varying the parameters as to obtain a Δχ 2 (C. L. = 68%, p = 4) = 4.7. The amplitude and period of the long-term modulation we found are similar to the values measured by Sanna et al. (2017) from an analysis of the outbursts observed until 2015. The large χ 2 of the fit suggests caution in interpreting these results. SAX J1808 orbital variability is similar to that observed in black widow and redback millisecond pulsars, rotation-powered pulsars in close binary orbits that ablate matter from their very low-mass (1 M e ) companion star. Yet the presence of a sinusoidal-like modulation of the orbital phase and of a much lower, formally negative, orbital period derivative evolution than previously estimated appear to be solid enough conclusions to draw. The sinusoidal modulation is hardly explained through the presence of a third body. The mass function (see, e.g., Bildsten & Chakrabarty 2001) of a putative third body would be ;2.7 × 10 −8 M e . Considering a NS mass of ∼1.4 M e and neglecting the companion mass (;0.05 M e ), the implied mass for the hypothetical third body would be ∼0.004 M e , for a third body inclination similar to the one of the system, i ∼ 69° (Goodwin et al. 2019). However, assuming that the orbit of SAX J1808 and of the putative third body are planar, the expected Doppler modulation of the pulsar frequency is δν ∼ (2π/P) A ν ∼ 42 μHz, which is about 2 orders of magnitude higher than observed (see Figure 2).

Discussion
The long-term orbital evolution of SAX J1808 has been discussed by several authors (see, e.g., Di Salvo et al. 2008;Hartman et al. 2008;Burderi et al. 2009;Hartman et al. 2009;Patruno et al. 2012Patruno et al. , 2017Sanna et al. 2017). A conservative mass transfer was soon excluded as the mass accretion rate implied by the large   P 4 10 s s orb 12 1 --, indicated by the first outbursts, is 2 orders of magnitude larger than ≈10 −12 M e yr −1 estimated from the average X-ray flux observed summing outbursts and quiescence periods (Marino et al. 2019). Di Salvo et al. (2008) and Burderi et al. (2009) discussed the surprisingly large value of  P orb of SAX J1808 in terms of mass lost by the system at a rate of ≈10 −9 M e yr −1 from the inner Lagrangian point, e.g., due to irradiation by a rotation-powered pulsar active in quiescence . As also noted by Hartman et al. (2008) and Di Salvo et al. (2008), the fast orbital evolution of SAX J1808 is reminiscent of black widow and redback pulsars. In these systems, the orbital period may change unpredictably with time, with T asc variations ranging from a few seconds to a few tens of seconds over a timescale of tens of years (see, e.g., Ridolfi et al. 2016;Freire et al. 2017;Kumari et al. 2022). The black widow PSR J2051−0827 exhibits orbital variability characterized by sinusoidal modulation with changing amplitude (see Figure 5 from Shaifullah et al. 2016). A chaotic orbital evolution has been also observed in the transitional redback system PSR J1023 + 0038 during its radio pulsar state (Archibald et al. 2015), while the orbital period variations seem to be smoother in the X-ray active state (Jaodand et al. 2016;Papitto et al. 2019;Burtovoi et al. 2020;Illiano et al. 2023). The long-term orbital modulation of a few black widow pulsars (Applegate & Shaham 1994;Arzoumanian et al. 1994;Doroshenko et al. 2001) has been interpreted in terms of gravitational quadrupole coupling (GQC) model (Applegate 1992;Applegate & Shaham 1994). This model was suggested to apply also to the case of SAX J1808 by Patruno et al. (2012; see also Patruno et al. 2017;Sanna et al. 2017). It envisages a gravitational coupling between the orbit and variations in the companion quadrupole moment, ΔQ, due to cyclic spin-up and spin-down of the star's outer layers. If ΔQ > 0, the companion will become more oblate, its gravitational potential in the equatorial plane will increase, and the orbit will shrink (  P 0 orb < ). On the contrary, if ΔQ < 0, the companion star will become less oblate, and the orbit will expand (  P 0 orb > ). Torques produced by magnetic activity of the companion would generate the angular momentum variations that are rapidly transmitted to the orbit. Blue points represent the measurements made with RXTE from the 1998 outburst to that of 2011 (Hartman et al. 2008;Burderi et al. 2009;Papitto et al. 2011), the green dot is the best value found for the 2015 outburst (Sanna et al. 2017), the orange one for the 2019 outburst (Bult et al. 2020), and the red one is from this work. The dotted line indicates a quadratic fitting function, while the dashed line indicates the best-fitting quadraticsinusoidal function. Middle panel: residuals relative to the quadratic model. Bottom panel: residuals relative to the quadraticsinusoidal fitting function. We point out that different y-axes are used for the second and the third panels.
Given the observed parameters of the long-term oscillation of SAX J1808, the GQC mechanism requires the companion to feature a magnetic field with a strength of B 2 ; 6 × 10 3 G and provide an internal luminosity of L GQC ; 10 30 erg s −1 (see Equations (15) and (16) in Sanna et al. 2017, derived from Applegate 1992Applegate & Shaham 1994), taking the NS and the companion masses equal to M NS ; 1.4 M e and M 2 ; 0.05 M e , respectively. However, identifying the energy source required to power such a mechanism in the case of SAX J1808 is not trivial, since nuclear burning or external irradiation of the companion star cannot provide such a high luminosity (Patruno et al. 2017;Sanna et al. 2017). Sanna et al. (2017) proposed that tidal dissipation could provide such a power if the donor is maintained in asynchronous rotation compared to the orbit by a magnetic braking mechanism. Irradiation by the pulsar wind would sustain the relatively high mass-loss rate required. In fact, in order to provide the required L GQC , the secondary would have to lose mass at a rate (Applegate & Shaham 1994;Sanna et al. 2017): = -yr −1 ), α = ℓ ej /ℓ 2 represents the amount of specific angular momentum that is carried away by such an outflow in units of the secondary specific angular momentum, β is the fraction of mass lost by the companion that is accreted onto the NS, and g(β, q, α) = 1 − β q − (1 − β) (α + q/ 3)/(1 + q).
First, considering a magnetic lever arm l ; 0.5a (similarly to Applegate & Shaham 1994;Sanna et al. 2017), the mass-loss rate is estimated to be   m 1.7 9 -(Equation (4)). Assuming that only a fraction β = 0.01 of the mass transferred by the companion is accreted by the NS, while the rest is ejected with the specific angular momentum at the inner Lagrangian point  (5)). Such a positive derivative seems too large to be compatible with the orbital phase evolution we found. Fixing the  P orb in Equation (3) to such a large value and repeating the fit leads to an unreasonably high value of the fit χ 2 (15817.9/4). Second, assuming a magnetic lever arm l = a (in analogy with what was done in Applegate & Shaham 1994), the mass-loss rate is estimated to be   m 0.4 9 -(Equation (4)). For α ; 0.7, we obtain  P 1.6 orb, 12 = -(Equation (5)), still too large for the observed orbital evolution (χ 2 /dof = 808.5/4).
By considering a range of orbital period derivative  P orb, 12 between 0 and −0.55 (i.e., ± four times the uncertainty on the best-fitting value) we deduce a range of values for α between 1.03 and 1.06 (see Equation (5)) for a mass-loss rate of   m 1.7 9 -(l ; 0.5a). If we assume   m 0.4 9 -(for l = a), the range of value for a is between 1.02 and 1.13.
While previous models had to assume that mass left the binary system with the specific angular momentum at the inner Lagrangian point (in order to yield a large positive orbital period derivative, see Equation (5)), the analysis of the data set presented here suggests that the orbit is contracting at a rate 1 order of magnitude lower than the expansion previously reported. As a consequence, the mass has to leave the system with an angular momentum equal to or greater than that of the secondary center of mass, so as to make the last term in Equation (5) smaller than the first term. A magnetic slingshot mechanism (see, e.g., Ferreira 2000;Waugh et al. 2021;Faller & Jardine 2022) by the strong B-field (B 2 ; 6 × 10 3 G) of the companion required to power the observed GQC luminosity might contribute to increase the specific angular momentum carried away by the matter ejected by the pulsar wind from the inner Lagrangian point. The observations of future outbursts will confirm the parameters of the long-term sinusoidal modulation, and help constrain the sign and magnitude of the orbital period derivative, which largely influence the conclusions on the rate of mass loss required to power the GQC mechanism and the location from which mass is ejected.

Conclusions
We presented a coherent timing analysis of NICER observations of SAX J1808.4−3658 during its 2022 outburst. We updated the orbital solution and investigated the pulse phase evolution during the outburst. We focused on the fundamental frequency, since the second harmonic was often too weak to be detected. We first modeled the phase delays using a constant frequency model, because the addition of a quadratic term (i.e.,  0 n ¹ ) did not produce a significant improvement in the data description. Because of the still large phase residuals, we then added to the linear model a dependence of the pulse phase on the flux, following Bult et al. (2020), significantly improving the fit's χ 2 . We observed an anticorrelation between the phase delays and the source flux, that holds only for count rates lower than ∼100 c/s, i.e., in the reflaring phase.
We confirmed the secular spin-down of   10 SD 15 n --Hz s −1 , as found in previous works (see, e.g., Patruno et al. 2012;Sanna et al. 2017;Bult et al. 2020), compatible with the energy losses expected from a ≈10 26 G cm 3 rotating magnetic dipole.
For the first time in the last twenty years, the orbital phase evolution showed evidence that the orbit has contracted since the last epoch. The long-term behavior of the orbit is described by a ∼11 s modulation with a ∼21 yr period. We excluded the presence of a third body, as the expected Doppler modulation