Timing Analysis of the 2022 Outburst of the Accreting Millisecond X-Ray Pulsar SAX J1808.4-3658: Hints of an Orbital Shrinking

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_⊙ 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³⁶ erg s⁻¹ (Gilfanov et al. 1998) at the peak of the outbursts, and decreases down to L_X ∼ few × 10³¹ erg s⁻¹ 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³¹ erg s⁻¹, and L_UV ≈ 2 × 10³² erg s⁻¹) 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 long-term spin and orbital evolution of the pulsar and discuss the implications for the models of AMXPs.

The NASA X-ray telescope Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2012) monitored SAX J1808 from 2022 August 19 (MJD 59810) until 2022 October 31 (MJD 59883) (ObsIDs starting with 505026 and 557401). The top panel of Figure 1 shows the 0.5–10 keV light curve. Visibility constraints prevented NICER from obtaining a homogeneous coverage of the outburst. The data were reduced and processed using HEASoft version 6.30 and nicerl2 task (NICERDAS version 7a), retaining events in the 0.5–10 keV energy range. We corrected the photon arrival times to the Solar System Barycenter (SSB) using the JPL ephemerides DE405 (Standish 1998). We adopted the source coordinates R.A. (J2000) = 18^h08^m27647(7) and decl. (J2000) = $-36^\circ 58^{\prime} 43\buildrel{\prime\prime}\over{.} 90(25)$ from Bult et al. (2020). We estimated the background contributions to our data with the nibackgen3C50 tool (Remillard et al. 2022).

Zoom In Zoom Out Reset image size

The NICER monitoring started when SAX J1808 had almost attained its peak count rate of ∼300 c s⁻¹ (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 (χ²/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³⁶ erg s⁻¹ (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 short-lived 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).

3.1. 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 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):

$$\begin{eqnarray}&&{\rm{\Delta }}\phi (t)={\phi }_{0}-{\rm{\Delta }}\nu \left(t-{T}_{0}\right)-\displaystyle \frac{1}{2}\,\dot{\nu }\,{\left(t-{T}_{0}\right)}^{2}+{R}_{\mathrm{orb}}(t).\end{eqnarray} \tag{ 1 }$$

Here, ν is the pulsar spin frequency, T₀ is the chosen reference epoch, ϕ₀ is the pulse phase at T₀, Δν = ν(T₀) − ν_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 χ² 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 ${\sigma }_{{\nu }_{\mathrm{pos}}}\leqslant {\nu }_{0}\,y\,{\sigma }_{\gamma }\,{(1+{\sin }^{2}\beta )}^{1/2}\,2\pi /{P}_{\oplus }$, 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 ${\sigma }_{{\nu }_{\mathrm{pos}}}\leqslant 5\times {10}^{-8}\,\mathrm{Hz}$. We added in quadrature this systematic uncertainty to the statistical error of the spin frequency reported in Table 1.

Table 1. Timing Solution for SAX J1808 2022 Outburst

Parameter	Value
Epoch (MJD)	59810.5956860
${a}_{1}\sin i$ (lt-s)	0.0628033(57)
Porb (s)	7249.1600(13)
Tasc (MJD)	59810.6179996(17)
Linear phase model
ν (Hz)	400.975209557(50)
χ2/dof	699.1/285
Flux-adjusted phase model
ν (Hz)	400.975209535(50)
b	1.44(49)
Γ	−0.81(12)
χ2/dof	450.0/283

Note. The timing solution was obtained adopting the source coordinates from Bult et al. (2020). Uncertainties are the 1σ statistical errors.

Download table as:  ASCIITypeset image

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 (Patruno et al. 2009; Bult et al. 2020). We modeled the phase delays using either a constant frequency model (i.e., setting $\dot{\nu }=0$ in Equation (1); see Table 1) or a constant spin frequency derivative. The quadratic fit returns a value of the average frequency derivative of $\dot{\nu }$ = 2.4(4.0) × 10⁻¹⁵ Hz s⁻¹ (χ²/dof = 698.2/284), which is compatible with zero. The probability of a chance improvement of the χ² 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 $\dot{\nu }=4.4(8)\times {10}^{-13}$ Hz s⁻¹ 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⁻¹³ Hz s⁻¹ for a 1.4 M_⊙ NS accreting at a rate of 10⁻⁹ M_⊙ yr⁻¹ from a disk truncated at 20 km from the NS; see, e.g., Di Salvo et al. 2019 and references therein). Hartman et al. (2008, 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 ($| \dot{\nu }| \lt 2.5\times {10}^{-14}$ Hz s⁻¹). 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, ${R}_{\mathrm{flux}}{(t)={{bF}}_{X}(t)}^{{\rm{\Gamma }}}$, where F_X is the X-ray bolometric flux, b = −0.87(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 ${\dot{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 R(t) (R_flux(t) = bR(t)^Γ). The resulting χ² shows a significant improvement with respect to the linear phase model (F-test probability of ∼8.5 × 10⁻²⁸; 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⁻¹³ Hz s⁻¹ (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⁻¹, 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 (δ ϕ ∼ 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.

3.2. 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 long-term 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 ${\dot{\nu }}_{\mathrm{SD}}\simeq -{10}^{-15}$ Hz s⁻¹, compatible with the energy losses expected from a ≈10²⁶ G cm³ 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⁻⁷ 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:

$$\begin{eqnarray}&&{\rm{\Delta }}\nu (t)=\delta {\nu }_{98}+{\dot{\nu }}_{\mathrm{SD}}(t-{T}_{98})+\delta {\nu }_{\mathrm{pos}}(t,\lambda ,\beta ).\end{eqnarray} \tag{ 2 }$$

Here, δ ν₉₈ = ν(t) − ν₉₈ is the spin frequency difference compared to the 1998 value, ν₉₈ = 400.975210371 Hz (Hartman et al. 2008), T₉₈ = 50914.8 MJD (Hartman et al. 2008), and δ ν_pos(t, λ, β) is the Doppler correction (see, e.g., Bult et al. 2020). We found δ ν₉₈ = 2.7(1.9) × 10⁻⁸ Hz, ${\dot{\nu }}_{\mathrm{SD}}=-1.152(56)\times {10}^{-15}$ Hz s⁻¹, δ λ = 042(15), and δ β = −093(38), with χ²/dof =34.9/5. Uncertainties of our best-fitting values were estimated from the parameters' range required to increase the χ² from the fit by an amount Δχ²(C. L. = 68%, p = 4) = 4.7, where p is the number of interesting free parameters (Avni 1976; Lampton et al. 1976; Yaqoob 1998). The spin-down trend observed across the previous outbursts is therefore confirmed. Also, the coordinate offsets are compatible within 1σ with what was found by Bult et al. (2020), and correspond to R.A.(J2000) = 18^h08^m27656(12), decl. (J2000) = $-36^\circ 58^{\prime} 44\buildrel{\prime\prime}\over{.} 222(89)$.

Zoom In Zoom Out Reset image size

3.3. Orbital Period Evolution

To investigate the orbital evolution, we computed the difference ΔT_asc between the measurements of the epoch of passage at the ascending node during the various outbursts and the values extrapolated from the epoch of passage at the ascending node estimated in the 2002 outburst (T_ref = 52499.9602472 MJD), assuming a constant orbital period (P_ref = 7249.156980(4) s; Hartman et al. 2009), ΔT_asc,i = T_asc,i − (T_ref + N_orb P_ref). Here, T_asc,i is the epoch of passage at the ascending node for the ith outburst, and N_orb is the nearest integer number of orbital cycles since T_ref. Until the 2008 outburst, the orbital phase evolution was consistent with an expansion at an average rate of ≃4 × 10⁻¹² s s⁻¹ (Di Salvo et al. 2008; Hartman et al. 2008; Burderi et al. 2009; Hartman et al. 2009). Subsequent outbursts first suggested an acceleration of the expansion (Patruno et al. 2012), then a transition to a slower evolution (Patruno et al. 2017; Sanna et al. 2017; Bult et al. 2020). The orbital phase we measured in 2022 data (see the red point in the top panel of Figure 3) indicates the first decrease of the orbital period seen from SAX J1808 in the last 20 yr. Indeed, modeling the ΔT_asc evolution with a constant period derivative (dotted line in Figure 3), leaves evident residuals with a sinusoidal shape (χ²/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:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{T}_{\mathrm{asc}}({N}_{\mathrm{orb}})=\delta {T}_{\mathrm{ref}}+\delta {P}_{\mathrm{ref}}{N}_{\mathrm{orb}}\\ \quad +\,\displaystyle \frac{1}{2}\,{\dot{P}}_{\mathrm{orb}}\,{P}_{\mathrm{ref}}\,{N}_{\mathrm{orb}}^{2}\\ \,+\,A\sin \left[\displaystyle \frac{2\,\pi }{P}\,({N}_{\mathrm{orb}}-N)\right].\end{array}\end{eqnarray} \tag{ 3 }$$

Here, δ T_ref is the offset from the 2002 epoch of passage at the ascending node, δ P_ref is the correction to the orbital period at the epoch of the 2002 outburst, ${\dot{P}}_{\mathrm{orb}}$ is the orbital period derivative, and A, P, and N are the amplitude, period, and phase of the sinusoidal function, respectively. The addition of the last term in Equation (3) led to a decrease of the fit's χ²/dof down to 117.9/3. Although statistically speaking the fit is clearly still unacceptable, an F-test indicates that the probability that the improvement occurs by chance is 0.1%. The best-fit values are the following: δ P_ref = 4.63(16) × 10⁻⁴ s, and ${\dot{P}}_{\mathrm{orb}}=-2.82(69)\,\times {10}^{-13}\,{\rm{s}}\ \,\ {{\rm{s}}}^{-1}$, A = 11.30(33) s, and P = 7.57(21) × 10³ days. We evaluated the uncertainties by varying the parameters as to obtain a Δχ²(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 χ² 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_⊙) 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⁻⁸ M_⊙. Considering a NS mass of ∼1.4 M_⊙ and neglecting the companion mass (≃0.05 M_⊙), the implied mass for the hypothetical third body would be ∼0.004 M_⊙, 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).

Zoom In Zoom Out Reset image size

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. 2012, 2017; Sanna et al. 2017). A conservative mass transfer was soon excluded as the mass accretion rate implied by the large ${\dot{P}}_{\mathrm{orb}}\simeq 4\times {10}^{-12}\,{\rm{s}}\ \,\ {{\rm{s}}}^{-1}$, indicated by the first outbursts, is 2 orders of magnitude larger than ≈10⁻¹² M_⊙ yr⁻¹ 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 ${\dot{P}}_{\mathrm{orb}}$ of SAX J1808 in terms of mass lost by the system at a rate of ≈10⁻⁹ M_⊙ yr⁻¹ from the inner Lagrangian point, e.g., due to irradiation by a rotation-powered pulsar active in quiescence (Burderi et al. 2003). 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 (${\dot{P}}_{\mathrm{orb}}\lt 0$). On the contrary, if ΔQ < 0, the companion star will become less oblate, and the orbit will expand (${\dot{P}}_{\mathrm{orb}}\gt 0$). Torques produced by magnetic activity of the companion would generate the angular momentum variations that are rapidly transmitted to the orbit.

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₂ ≃ 6 × 10³ G and provide an internal luminosity of L_GQC ≃ 10³⁰ erg s⁻¹ (see Equations (15) and (16) in Sanna et al. 2017, derived from Applegate 1992; Applegate & Shaham 1994), taking the NS and the companion masses equal to M_NS ≃ 1.4 M_⊙ and M₂ ≃ 0.05 M_⊙, 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):

$$\begin{eqnarray}&&{\dot{m}}_{2}=0.25\times {10}^{-9}\,{\left(\displaystyle \frac{a}{l}\right)}^{2}\,{L}_{30}^{1/2}\,{t}_{\mathrm{sync},4}^{-1/2}\,{I}_{2,51}^{1/2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 4 }$$

Here, a is the orbital separation, l is the magnetic lever arm of the mass ejected from the companion star, L₃₀ is the tidal luminosity in units of 10³⁰ erg s⁻¹, I_2,51 is the companion moment of inertia in units of 10⁵¹ g cm², and t_sync,4 is the tidal synchronization time in units of 10⁴ yr. For a Roche lobe filling companion, the latter is estimated as ${t}_{\mathrm{sync}}=0.65\times {10}^{4}{{\mu }_{12}}^{-1}{(1+q)}^{2}\,{M}_{2,\odot }{{R}_{2,\odot }}^{-1}$ yr (Applegate & Shaham 1994), where ${\mu }_{12}=3\,{L}_{2,\odot }^{1/3}\,{R}_{2,\odot }^{-5/3}\,{M}_{2,\odot }^{2/3}$ is the mean dynamic viscosity in units of 10¹² g cm⁻¹ s⁻¹, and L_2,⊙ = L_GQC/L_⊙, M_2,⊙ and R_2,⊙ are the luminosity, the mass, and the radius of the companion star in Solar units, respectively. Assuming M₂ ≈ 0.05 M_⊙, R₂ ≈ 0.13 R_⊙ (Bildsten & Chakrabarty 2001), we obtain t_sync ≃ 3.4 × 10³ yr, similar to the value reported by Sanna et al. (2017). The corresponding variation of the orbital period in units of 10⁻¹² s s⁻¹ is expressed by Di Salvo et al. (2008); Burderi et al. (2009); Sanna et al. (2017):

$$\begin{eqnarray}\begin{array}{rcl}{\dot{P}}_{\mathrm{orb},-12} & = & -1.38\,{M}_{\mathrm{NS},\odot }^{5/3}\,q\,{(1+q)}^{-1/3}\,{P}_{\mathrm{orb},2h}^{-5/3}\\ & & +0.648\,{M}_{\mathrm{NS},\odot }^{-1}\,{q}^{-1}\,{P}_{\mathrm{orb},2h}\,g(\beta ,q,\alpha )\,{\dot{m}}_{-9}.\end{array}\end{eqnarray} \tag{ 5 }$$

Here, P_orb,2h is the orbital period in units of 2 h, ${\dot{m}}_{-9}={\dot{m}}_{2}/({10}^{-9}\,{M}_{\odot }$ yr ⁻¹), α = ℓ_ej/ℓ₂ 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 ${\dot{m}}_{-9}\simeq 1.7$ (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 ($\alpha ={[1-0.462\,{q}^{1/3}\,{(1+q)}^{2/3}]}^{2}\simeq 0.7$; Di Salvo et al. 2008; Burderi et al. 2009; Sanna et al. 2017) requires a period derivative of ${\dot{P}}_{\mathrm{orb},-12}=7.0$ (Equation (5)). Such a positive derivative seems too large to be compatible with the orbital phase evolution we found. Fixing the ${\dot{P}}_{\mathrm{orb}}$ in Equation (3) to such a large value and repeating the fit leads to an unreasonably high value of the fit χ² (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 ${\dot{m}}_{-9}\simeq 0.4$ (Equation (4)). For α ≃ 0.7, we obtain ${\dot{P}}_{\mathrm{orb},-12}=1.6$ (Equation (5)), still too large for the observed orbital evolution (χ²/dof = 808.5/4).

By considering a range of orbital period derivative ${\dot{P}}_{\mathrm{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 ${\dot{m}}_{-9}\simeq 1.7$ (l ≃ 0.5a). If we assume ${\dot{m}}_{-9}\simeq 0.4$ (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₂ ≃ 6 × 10³ 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.

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., $\dot{\nu }\ne 0$) 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 χ². 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 ${\dot{\nu }}_{\mathrm{SD}}\simeq -{10}^{-15}$ Hz s⁻¹, 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²⁶ G cm³ 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 of the pulsar frequency would be about 2 orders of magnitude higher than observed.

We discussed the observed ${\dot{P}}_{\mathrm{orb}}=-2.82(69)\times {10}^{-13}\,{\rm{s}}\ \,\ {{\rm{s}}}^{-1}$ in terms of a coupling between the orbit and variations in the mass quadrupole of the companion star (GQC model; Applegate 1992; Applegate & Shaham 1994). Data suggest that matter leaving the system with the specific angular momentum of the companion center of mass could maintain the donor star out of tidal locking and drive the required oscillation of its structure. A strong magnetization of the companion star (B ≃ 6 × 10³ G at the surface) is required to couple the mass loss to the donor star's rotation and to increase the angular momentum carried away by the ejected matter compared to the orbital value.

Based on past recurrence times, it is expected that there will be a new outburst of SAX J1808 in approximately 3 yr (2025). The observations of the next outburst will be of paramount importance to confirm the source's orbital evolution, by decreasing the correlation between the long-term modulation of the orbital phase epoch and the quadratic term that represents a secular orbital period derivative. This would constrain even more the mass-loss rate and the location from which mass is ejected needed to power the GQC mechanism. Detecting pulsations during the quiescent state would greatly increase our ability to track the pulsar orbital evolution without relying solely on data taken during unpredictable outbursts. Even though a rotation-powered pulsar is expected to turn on during quiescence (Burderi et al. 2003), deep searches for radio (Burgay et al. 2003; Patruno et al. 2017) or gamma-ray (de Oña Wilhelmi et al. 2016) pulsations have not succeeded in detecting a signal, so far. Recently, the discovery of optical pulsations from a couple of millisecond pulsars (Ambrosino et al. 2017, 2021; A. Miraval Zanon et al. 2022, in preparation) opened the intriguing possibility of exploiting the higher sensitivity in this band compared to higher energies, and we plan to use this additional diagnostic also to investigate the orbital evolution of this source.

This work is based on observations acquired with the NASA mission NICER. This research has made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. G.I. is supported by the AASS Ph.D. joint research program between the University of Rome "Sapienza" and the University of Rome "Tor Vergata" with the collaboration of the National Institute of Astrophysics (INAF). F.A., G.I., A.M.Z., A.P., L.S., and D.d.M. acknowledge financial support from the Italian Space Agency (ASI) and National Institute for Astrophysics (INAF) under agreements ASI-INAF I/037/12/0 and ASI-INAF n.2017-14-H.0, from INAF Research Grant "Uncovering the optical beat of the fastest magnetised neutron stars (FANS)". F.A., G.I., A.M.Z., A.P., and L.S. also acknowledge funding from the Italian Ministry of University and Research (MUR), PRIN 2020 (prot. 2020BRP57Z) "Gravitational and Electromagnetic-wave Sources in the Universe with current and next-generation detectors (GEMS)". L.S. and T.D.S. acknowledge financial contributions from "iPeska" research grant (P.I. Andrea Possenti) funded under the INAF call PRIN-SKA/CTA (resolution 70/2016). L.S. is also partially supported by PRIN-INAF 2019 no. 15. T.D.S. acknowledges financial support from PRIN-INAF 2019 (n. 89, PI: Belloni). A.M.Z. is supported by PRIN-MIUR 2017 UnIAM (Unifying Isolated and Accreting Magnetars, PI S. Mereghetti). P.B. acknowledges support from NASA through CRESST II cooperative agreement (80GSFC21M0002). F.C.Z. and A.M. are supported by the H2020 ERC Consolidator Grant "MAGNESIA" under grant agreement No. 817661 (PI: Rea) and National Spanish grant PGC2018-095512-BI00. This work was also partially supported by the program Unidad de Excelencia Maria de Maeztu CEX2020-001058-M, and by the PHAROS COST Action (No. CA16214). F.C.Z. is also supported by Juan de la Cierva Fellowship. J.P. acknowledges support from the Academy of Finland grant 333112.

Facilities: ADS - , HEASARC - , NICER - .

Software: HEASoft (v6.30), NICERDAS (v7a), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020).