© The Authors 2026

1. Introduction

Accreting millisecond pulsars (AMSPs) are weakly magnetized (B ∼ 10^8 − 9 G) neutron stars (NSs) that accrete matter from a low-mass companion (≲1 M_⊙) via an accretion disk (e.g., Patruno & Watts 2021; Di Salvo & Sanna 2022; Campana & Di Salvo 2018). Their high spin frequencies result from a previous billion-year-long evolutionary phase of mass accretion (Alpar et al. 1982; Bhattacharya & van den Heuvel 1991). In this scenario, known as “recycling”, slowly rotating radio pulsars are spun up to millisecond periods through the accretion of matter and angular momentum transferred from the companion star via Roche-lobe overflow. AMSPs are all transient systems that spend most of their life in quiescence. Occasionally, they start several-week-long mass-accretion outbursts, during which the X-ray luminosity increases by a few orders of magnitude up to ≃10^{36 − 37} erg s⁻¹. During these outbursts, the pulsar magnetosphere channels the accreting matter towards the NS magnetic poles, producing coherent X-ray pulsations at a frequency higher than ∼100 Hz. More than two dozen AMSPs have been discovered so far (see, e.g., the most recent discovery in Ng et al. 2024).

The radius at which the NS magnetosphere starts controlling the dynamics of the infalling matter, R_in, decreases with the mass-accretion rate (Ghosh & Lamb 1979; Spruit & Taam 1993). Toward the end of an X-ray outburst, R_in expands and might exceed the corotation radius, R_co. When this happens, a centrifugal barrier is expected to hinder mass accretion onto the NS, in the so-called propeller regime (Illarionov & Sunyaev 1975; Stella et al. 1994; Campana et al. 1998). However, Spruit & Taam (1993) showed that accretion may proceed at a reduced rate even when R_in ≳ R_co, because the centrifugal barrier is too weak to expel matter from the system. In this regime, a “trapped” or “dead” disk may form (D’Angelo & Spruit 2010; D’Angelo & Spruit 2012), where matter accumulates due to inefficient angular-momentum loss, allowing only a small fraction to reach the NS surface. Magnetohydrodynamic simulations by Romanova et al. (2018) showed that either a weak or strong propeller state may be realized. In the former (R_in ≳ R_co), part of the inflowing material can penetrate the centrifugal barrier and accrete along the magnetic-field lines. In the strong propeller regime (R_in ≫ R_co), most of the matter is expelled from the system, and accretion onto the stellar surface is largely suppressed.

As they fade into quiescence, X-ray observations of AMSPs are a powerful tool with which to test these scenarios. SAX J1808.4−3658 (hereafter SAX J1808) is the best studied AMSP during the 11 approximately one-month-long outbursts shown from its discovery in 1998 (Wijnands & van der Klis 1998). After reaching a peak luminosity of ∼10^{36 − 37} erg s⁻¹ a few days after the outburst onset, this source typically shows a slow decay (∼10 days) followed by a rapid drop (3–5 days) (see, e.g., Fig. 1). Before returning to a quiescent level of ∼5 × 10³¹ erg s⁻¹ (Campana et al. 2004; Stella et al. 2000), SAX J1808 displays peculiar luminosity variations dubbed “reflares” for a month or longer. During reflares, the X-ray luminosity varies in the range of ∼10^{32 − 35} erg s⁻¹ on timescales of approximately one-to-two days (Patruno et al. 2009a, 2016; Patruno & Watts 2021). It has been proposed that reflares originate from small changes in the disk density at the end of an outburst that increase the outer-disk temperature, partially ionizing hydrogen (Patruno et al. 2009a, 2016). This leads to a rapid increase in accretion into the inner disk and a subsequent temporary rebrightening. The detection of X-ray pulsations, the aperiodic variability and the stable spectral distribution of the X-ray emission observed during the reflares of SAX J1808 suggest that the innermost regions of the accretion flow remain close to the corotation radius even at the lowest X-ray luminosity attained (a few ×10³⁴ erg s⁻¹; Patruno et al. 2009a, 2016). Coherent X-ray pulsations have also been observed at similarly low luminosities of a few ×10³³ erg s⁻¹ in two additional AMSPs (Bult et al. 2019a; Illiano et al. 2026). Altogether, these observed properties suggest that channeled accretion in AMSPs can persist down to very low mass-accretion rates at which propeller ejection of matter would be rather expected.

Fig. 1. Light curve of 2022 outburst of SAX J1808 in 0.5–10 keV energy band, using 1 ks bins. NICER observations are shown in black, and XMM-Newton observations are given in red. The blue star marks the epoch of the HST/STIS observation (2022 September 10). We rescaled the XMM-Newton count rate to the NICER count rate using a conversion factor of 1.44, which was calculated with the WebPIMMS tool (https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl), and a power-law model with a photon index of Γ = 2.04 (Table 2).

SAX J1808 was also the first AMSP to show coherent pulsations at optical and UV wavelengths during the rising and final stages of its 2019 outburst, respectively (Ambrosino et al. 2021). The observed highly pulsed optical and UV luminosities ($L_{\mathrm{pulsed}}^{\mathrm{opt}}\approx 3\times {10}^{31}$ erg s⁻¹, $L_{\mathrm{pulsed}}^{\mathrm{UV}}\approx 2\times {10}^{32}$ erg s⁻¹) exceeded the values expected in the case of cyclotron emission from accretion columns by roughly 50 − 100 times, and the expected thermal output at those wavelengths by several orders of magnitude. Optical and UV pulsations have also been detected from the transitional millisecond pulsar PSR J1023+0038 (Ambrosino et al. 2017; Papitto et al. 2019; Jaodand et al. 2021; Ambrosino et al. 2021; Miraval Zanon et al. 2022), and from the redback PSR J2339-0533 (Papitto et al. 2025). In the former case, the high optically pulsed luminosities observed required a nonstandard interpretation in terms of synchrotron emission from the interface between the relativistic wind of a rotation-powered pulsar and the accretion flow (Papitto et al. 2019; Veledina et al. 2019; Baglio et al. 2023, 2025). However, this scenario assumes that the NS is active, as a rotation-powered pulsar and cannot apply to SAX J1808, where simultaneous X-ray pulsations indicate that the accretion flow extends down to the NS surface. The origin of the pulsed optical and UV emission in SAX J1808 thus remains uncertain.

In this paper, we report on a 125 ks-long XMM-Newton observation and a simultaneous 2.2 ks Hubble Space Telescope (HST) observation, performed during the reflaring stage of the outburst observed in SAX J1808 in 2022 (Illiano et al. 2023b). XMM-Newton’s uninterrupted coverage enabled a high-resolution study of the spectral and temporal evolution of the source X-ray emission as its flux varies by roughly one order of magnitude down to a few 10³⁴ erg s⁻¹. The simultaneous HST observation allowed us to extend the range of mass-accretion rates at which UV pulsations have been studied by a factor of approximately two and to perform the first simultaneous study in X-ray and UV bands.

2. Observations

2.1. XMM-Newton

XMM-Newton (Jansen et al. 2001) observed SAX J1808 for 125 ks starting on 2022 September 9 at 14:22:09 (UTC) (Obs.ID. 0884700801, PI: A. Papitto), during the final stage of the outburst (see red points in Fig. 1). Figure 2 shows the background-subtracted 0.5–10 keV EPIC-pn light curve. We used the Science Analysis System (SAS; v20.0.0) to process and reduce the data. We barycentered the photon arrival times observed by XMM-Newton using the barycen tool, with the JPL DE-405 Solar System ephemerides and the source position derived by Bult et al. (2020). The EPIC-pn camera operated in timing mode to achieve the 29.5 μs temporal resolution needed to study source variability and minimize photon pile-up, and it was equipped with a thick filter. The effective exposure was reduced to 108.2 ks after excluding soft-proton flaring episodes, identified when the EPIC-pn 10–12 keV count rate exceeded 0.8 c s⁻¹. In timing mode, spatial information along one of the optical axis was discarded to allow a faster readout. As a result, the source appears on the EPIC-pn detector as a bright strip rather than its usual circular shape. The maximum number of counts were recorded in pixels characterized by RAWX coordinates 37 and 38. EPIC-pn events were extracted from a 21-pixel-wide region (1 pixel ≃4.1′) around the source position, spanning from RAWX = 28 to 48. We estimated the background far from the source, in a three-pixel-wide region centered on RAWX = 4.

Fig. 2. Background-subtracted 0.5–10 keV XMM-Newton/EPIC-pn light curve of 2022 outburst of SAX J1808 measured using 200 s-long bins. A red band indicates the low-flux interval, and a blue band shows the interval of the HST observation discussed in this paper (see Sect. 3.2 for details). Two Type I X-ray bursts, detected at 59831.87 and 59832.82 MJD, are not plotted. Green dash-dotted lines indicate the epochs of their occurrence.

During the XMM-Newton pointing, SAX J1808 exhibited two Type I X-ray bursts at MJD 59831.87 and 59832.82, which will be discussed in a dedicated upcoming publication. After obtaining the light curve, we created good time intervals (GTIs), discarding time intervals starting 10 s before and ending 150 s after each burst onset to perform spectral and timing analyses of the non-bursting emission. Following the recommended XMM-Newton/SAS procedures¹, we extracted the source and background spectra from the cleaned event file, considering only single- and double-pixel events (PATTERN ≤ 4) and filtering out spurious events (FLAG = 0). We generated spectral response files using the tasks rmfgen and arfgen. The EPIC-pn spectra were rebinned using the SAS specgroup tool to have at least 25 background-subtracted counts per spectral channel and to sample the instrument’s resolution full width at half maximum with no more than three channels at any energy.

2.2. HST/STIS

The Space Telescope Imaging Spectrograph (STIS) on board HST observed SAX J1808 for 2.2 ks on 2022 September 10, starting at 02:05:41 (UTC) (GO/DD-17245, PI Miraval Zanon). The UV observation occurred during the interval of minimum X-ray flux, as shown by the blue band in Fig. 2. We removed an ∼400 s bump in the recorded flux at MJD 59832.103, which occurred when the observatory grazed the region of the South Atlantic Anomaly (HST help desk, private communication). The resulting UV light curve of the source observed with STIS is consistent with a constant count rate of ∼45 c s⁻¹. The spectroscopic observation was performed in TIME-TAG mode, achieving a time resolution of 125 μs. The G230L grating, equipped with a 52 × 0.2 arcsec² slit, provided a spectral resolution of ∼500 over the nominal (first-order) range. The total photon count rate was R_HST = (46.62 ± 0.14) c s⁻¹, with approximately 35% attributed to background radiation (BKG_HST = (16.66 ± 0.09) c s⁻¹).

We extracted the source photons across 19 slit channels (990 − 1010) and in the 165 − 310 nm wavelength range. This selection isolates the source signal, minimizes background contributions, and avoids noise arising from the G230L grating’s reduced response at the edges of the wavelength range. We estimated the background signal by selecting photons in the STIS slit channels outside the source region (200 − 900 and 1100 − 2000) in order to have high statistics and avoid the source contribution. We then averaged the resulting value and normalized it to the number of source slit channels. To analyze HST/STIS data, we first corrected the position of the slit channels using the Python custom-built external function stis_photons². UV photon arrival times were then corrected to the Solar System barycenter (SSB) using the ODELAYTIME task (subroutine available in the IRAF/STDAS software package) and JPL DE200 ephemerides. Note that the difference between the JPL DE200 and DE405 ephemerides is negligible compared to the ∼1 s uncertainty on the absolute timing of the HST data (HST help desk, private communication).

2.3. NICER

We analyzed the Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2012) observations obtained during the 2019 and 2022 outbursts and previously presented by Bult et al. (2020) and Illiano et al. (2023b). In 2019, NICER observed SAX J1808 from July 30 (MJD 58695) until November 8 (MJD 58795; ObsIDs 205026 and 258401). In 2022, we observed the source between August 19 (MJD 59810) and October 31 (MJD 59883; ObsIs 505026 and 557401). We processed and corrected the data according to the procedures described by Bult et al. (2020) and Illiano et al. (2023b). All data were processed using the standard NICER analysis software NICERDAS, part of HEASoft (v6.35.2). We processed the data using the nicerl2 task, selecting events in the 0.5 − 10 keV energy range. Photon arrival times for both datasets were corrected to the Solar System barycenter using the JPL DE405 ephemerides and the BARYCORR tool. We then adopted the source coordinates from Hartman et al. (2008) and Bult et al. (2020), respectively, for the 2019 and 2022 NICER datasets, to allow a direct comparison with the analyses presented by Bult et al. (2020) and Illiano et al. (2023b).

3. Coherent timing analysis

3.1. X-ray timing analysis

We performed a coherent timing analysis of the 2022 outburst of SAX J1808 using the XMM-Newton/EPIC-pn high-time-resolution dataset. We first corrected the photon arrival times for the Doppler delays introduced by the binary motion using the orbital parameters measured by Illiano et al. (2023b) from the pulse-timing analysis of 2022 NICER data. We created pulse profiles by epoch folding 1000 s time intervals in 16 phase bins around the best available estimate of the spin frequency (ν_F = 400.975209557 Hz from Illiano et al. 2023b). We adopted a bin length of 1000 s to ensure adequate statistics to detect the pulsation and track its amplitude and phase variations (see Fig. 4). Two harmonic components at most are usually required to satisfactorily model the pulse profile of SAX J1808 during such short intervals (see, e.g., Hartman et al. 2008), justifying the choice of sampling it in 16 phase bins. Only statistically significant pulse profiles were selected, namely folded profiles with a ratio between the sinusoidal amplitude and the corresponding 1σ error equal to or greater than three. The second panel of Fig. 4 shows the fractional amplitudes of the first and second harmonic. We focused the analysis of this work on the former, as it is significantly stronger. During the minimum of the X-ray flux (MJD 59832.02 to 59832.14; see the red band in Fig. 2), the amplitude of the fundamental (black dots) roughly doubled up to a value of ∼10%. Figure 3 shows the X-ray pulse profiles obtained by epoch by folding the high and low-flux intervals (red band in Fig. 2) of the XMM-Newton observation separately.

Fig. 3. X-ray pulse profiles obtained by folding the high- and low-flux intervals of XMM-Newton observation into 16 phase bins at a spin frequency of νF = 400.975209557 Hz. The solid line represents the best-fit model (constant plus two harmonics), while the dashed red and blue lines indicate the first and second harmonic contributions, respectively. For clarity, two cycles are shown.

Fig. 4. Top panel: 0.5–10 keV XMM-Newton/EPIC-pn light curve using 1 ks bins, with Type I X-ray bursts removed. The dashed blue line shows the spline interpolation obtained with the UnivariateSpline function from scipy.interpolate. The blue band indicates the interval of the HST observation discussed in this paper (see Sect. 3.2 for details). Second panel: Pulse fractional amplitude for the first harmonic (black dots) and the second harmonic (red dots). Third and fourth panels: Phase residuals relative to a linear model and a model that includes a phase-flux correlation term, respectively. Note that, where not visible, the error bars are smaller than the data points.

We first modeled the evolution of the pulse phase measured over the first harmonic using the following relation (e.g., Burderi et al. 2007; Sanna et al. 2022b; Papitto et al. 2007):

$$\begin{array}{c}\hfill \varphi (t)={\varphi }_0-\mathrm{\Delta }\nu (t-T_0)-\frac{1}{2}\dot{\nu }{(t-T_0)}^2+R_{\mathrm{orb}}(t).\end{array}$$(1)

Here, T₀ is the reference epoch for the timing solution, ϕ₀ is the pulse phase at T₀, Δν = ν(T₀)−ν_F represents the difference between the frequency at the reference epoch and the spin frequency used to epoch-fold the data, and $\dot{\nu }$ is the average spin-frequency derivative. The term R_orb models the phase residuals arising from differences between the adopted orbital parameters and the actual ones (see, e.g., Deeter et al. 1981). We first considered a constant frequency model ($\dot{\nu }=0$), obtaining a value of ν(T₀) that was not compatible with that expected according to the timing solution reached by Illiano et al. (2023b) (see Table 1).

Adding a spin frequency derivative term to the model reduced the residuals (Δχ² = 339 for one degree of freedom less), but the modeling remained statistically unsatisfactory (χ_r² = 11.86 over 91 degrees of freedom). The residuals with respect to either a linear or a quadratic model show large deviations (see the third panel of Fig. 4). This phenomenon often affects the evolution of the pulse phase of AMSPs and is commonly referred to as timing noise (Burderi et al. 2006; Papitto et al. 2007; Hartman et al. 2008, 2009; Patruno et al. 2009b). Timing noise is characterized by phase variations with amplitudes of up to several tenths of a cycle, affecting the measurement of the pulse frequency and its derivative. Here, a large phase jump of ∼0.4 cycles occurring around 59832 MJD, on a timescale of a few hours, was the most evident feature (see Fig. 4). This jump coincides with a decrease in the X-ray flux of the source by a factor of approximately three. When the flux returned to the initial values, the pulse phase also reverted to the estimates obtained before the jump.

To account for the phase shifts, we added a term describing a correlation between the X-ray flux and the pulse phase to Eq. (1), R_flux(t), which is possibly related to hot spot drifts (Patruno et al. 2009b):

$$\begin{array}{c}\hfill \varphi (t)={\varphi }_0-\mathrm{\Delta }\nu (t-T_0)+R_{\mathrm{orb}}(t)+R_{\mathrm{flux}}(t).\end{array}$$(2)

Following Bult et al. (2020), we used R_flux(t) = b (F_X/F₀)^Λ, with F_X being the X-ray bolometric flux, approximated by the 0.5–10 keV count rate for simplicity, and F₀ the initial flux value at t = T₀. To apply such a model, we interpolated the light-curve data using the UnivariateSpline function (available in Python’s scipy.interpolate library), as done by Illiano et al. (2023b). This method allows the flux values to be interpolated at the exact time points where phase data are available. The estimates of the X-ray flux so obtained (plotted with a dashed blue line in the top panel of Fig. 4) were used in Eq. (2) to fit the pulse phases, obtaining the residuals shown in the bottom panel of Fig. 4. By adding this flux-phase correlation term, we obtained a reduced χ² of 1.82 for 90 degrees of freedom (d.o.f.), showing a highly significant improvement over the quadratic phase model (Δχ² = 915 for the addition of one free parameter), with a probability of ∼1.5 × 10⁻³⁸ of it being due to chance, according to an F-test. The best-fitting parameters of the flux-phase correlation term are $b=0.9_{-0.4}^{+1.6}$ and Λ = −0.17(9). Even though the χ² value obtained is still formally statistically unacceptable, the phase residuals no longer showed the phase jump, nor did they indicate the presence of significant structures or trends (Fig. 4, bottom panel). Unlike the linear and quadratic models, the value of the spin frequency obtained from the flux-adjusted model (Table 1) was consistent within the errors with that expected according to the timing solution reported in Illiano et al. (2023b), thus supporting the validity of the corrected model.

To estimate the uncertainties on the parameters b and Λ, we considered the contour levels of the fit χ², χ_min² + 2.3, and χ_min² + 4.61 (Fig. 5), corresponding to confidence levels of 68% and 90%, respectively, for a fit with two free parameters of interest (Lampton et al. 1976; Avni 1976; Yaqoob 1998). The shape of the contour plots strongly suggests a correlation b × Λ ≃ const between the fit parameters.

Fig. 5. Contour levels of χ2 obtained by varying both parameters b and Λ in the flux-adjusted phase fit. The contour levels are shown for χmin2 + 2.3 (blue) and χmin2 + 4.61 (red), corresponding to confidence levels of 68% and 90%, respectively, for a fit with two parameters (Lampton et al. 1976). The black dot indicates the values of b and Λ corresponding to the minimum χ2, while the dashed black line represents the contour corresponding to the constant value given by the product b × Λ.

For comparison, we applied the same coherent timing technique to investigate the phase evolution of NICER data collected during the 2019 and 2022 outbursts. To avoid introducing significant uncertainties into the spline function, we filtered out data from the end of the outburst when the temporal gaps between consecutive observations exceeded three days. In both datasets, a flux-adjusted phase model again provided a better description of the observed phases than a linear or quadratic modeling. We obtained $b=0.{75}_{-0.14}^{+0.30}$ and Λ = −0.14(5) for the 2019 dataset, and $b=0.{08}_{-0.03}^{+0.07}$ and Λ = −0.37(14) for the 2022 observation. The values of Λ are consistent, within uncertainties, with the value derived from the XMM-Newton timing analysis. These results are consistent with those of Bult et al. (2020), where Λ was fixed at −1/5. In contrast, the discrepancy with the values reported by Illiano et al. (2023b) arises from the different time interval considered for the phase fitting.

Figure 6 shows the fractional pulse amplitudes in seven energy bands, computed in the low-flux phase of the XMM-Newton light curve (top panel), and in the high-flux phase of the X-ray emission (bottom panel).

Fig. 6. Fractional amplitude of fundamental harmonic of X-ray pulse profiles in seven energy bands (0.5–1, 1–2, 2–3, 3–4, 4–5, 5–7, 7–10 keV) for the low-flux part of the light curve (top panel) and for the high-flux emission (bottom panel).

The morphology we observed closely matches interval 4 of Bult et al. (2020) (see their Fig. 2), corresponding to the reflare phase of the 2019 outburst, which shows the highest amplitudes and a fractional amplitude that increases with energy before returning to near-initial values at ∼10 keV. The amplitudes we observed are generally higher than those reported by Bult et al. (2020), although the relative increase in pulsed fraction from ∼1 keV to 2–3 keV is consistent across the two outbursts, both showing an increase of approximately 40% relative to the initial values.

To compare the properties of the X-ray and UV pulsations, we measured the X-ray pulse amplitude in an interval strictly simultaneous to the HST observations (marked by a blue band in Fig. 4). In that interval, we obtained A_X = (9.9 ± 0.8)%. The 0.5–10 keV flux measured in that time interval using a spectral model described in Table 2 was $F_{0.5-10}=(4.{24}_{-0.10}^{+0.14})\times {10}^{-11}{\mathrm{erg}\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. The pulsed 0.5–10 keV X-ray luminosity was then evaluated as $L_{\mathrm{pulsed}}^{\mathrm{X}}=(A_{\mathrm{X}}/\sqrt{2})(4\pi d^2{F}_{0.5-10})=(4.4\pm 0.4)\times {10}^{33}{d}_{3.5}^2{\mathrm{erg}\mathrm{s}}^{-1}$, where d_3.5 is the distance to the source in units of 3.5 kpc (Galloway & Cumming 2006).

3.2. UV timing analysis

To measure the spin period observed during the 2.2 ks-long HST/STIS observation, we performed an epoch folding search using m = 8 phase bins and a period resolution of δ_PHST, EFS = P²/2mT_obs = 1.7 × 10⁻¹⁰ s, where T_obs is the total exposure time of the HST observation. We thus obtained a best-fitting period of P_HST, EFS = 2.4939194(6)×10⁻³ s. Following Leahy (1987), the 1σ uncertainty reported was calculated using the following equation: σ_P = (P²/2 T_obs) ⋅ 0.71(χ_max²/m − 1)^−0.63, where the maximum of the χ² value was found to be χ_max² ≃ 26.6. The period inferred from UV data is consistent, within uncertainties, with that obtained from X-ray data. We then folded the HST data at the best spin period. A single sinusoidal component describes the resulting UV pulsed profiles well, and an additional harmonic is not statistically justified, as it does not significantly improve the fit chi-square (see Fig. 7). To take the background into account, we evaluated a correction factor as N_γ(tot)/(N_γ(tot) − N_γ(bkg))≃1.556(5), where N_γ(tot) is the total number of detected photons and N_γ(bkg) the number of background photons. After applying this correction, the root-mean-square (r.m.s.) amplitude resulted as $A_{\mathrm{UV}}^{\mathrm{rms}}=(1.9\pm 0.5)\%$.

Fig. 7. Pulse profile obtained by folding the whole HST/STIS observation (GO/DD-17245, PI Miraval Zanon) into eight phase bins at a spin frequency of νF = 400.975209557 Hz. The solid line represents the best-fit model, consisting of a single sinusoidal component. For clarity, two cycles are shown.

4. X-ray spectral analysis

We analyzed the X-ray spectral evolution of SAX J1808 during the outburst’s flaring stage by measuring the hardness ratio (HR), which is defined as the ratio of the count rates observed in the (2–10 keV) and (0.7–2 keV) energy bands. Panels (a) and (b) of Fig. 8 show the variation of the HR over time and intensity, respectively. The source emission significantly softened during the phase of minimum X-ray count rate. The continuous decrease in the HR reached a minimum when the source flux was at its lowest. The HR then returned to its initial value, following the trend of the light curve.

Fig. 8. Panel (a): Evolution of HR of SAX J1808 during XMM-Newton observation of reflare phase of 2022 outburst. The HR was estimated from the (2–10 keV)/(0.7–2 keV) energy band. For comparison, the XMM-Newton light curve is shown in blue in the background. Panel (b): Hardness-intensity diagram (HID) of SAX J1808 obtained using XMM-Newton EPIC-pn observation and 200 s bins. The intensity is given in counts per second. The points corresponding to the minimum of the source X-ray flux on the light curve, between 59832.02 MJD and 59832.14 MJD, are marked in red.

To investigate the softening of the emission, we extracted spectra considering both the low-flux phase of the light curve (between 59832.02 MJD and 59832.14 MJD; see the red band of Fig. 2 and corresponding red data points in Fig. 8) and the high-flux emission, obtained by excluding the entire low-flux interval and its rise and decay phases; i.e., from MJD 59831.88 to 59832.23. We performed the spectral analysis using the X-ray spectral fitting package HEASARC XSPEC (Arnaud 1996) version 12.14.0. We used TBabs to model the interstellar photoelectric absorption, with the chemical abundances of Wilms et al. (2000) and the cross-sections reported in Verner et al. (1996). To avoid the impact of calibration uncertainties at low energies, we excluded data below 1 keV from the EPIC-pn spectral analysis³. Following previous works, we fixed the equivalent hydrogen column density at 0.21 × 10²² cm⁻² (Papitto et al. 2009; Di Salvo et al. 2019). We first modeled the 1–10 keV spectrum observed in the high-flux interval with a blackbody component, bbodyrad, and a thermal Comptonization component, nthComp (Zdziarski et al. 1996; Zycki et al. 1999). This spectral decomposition is commonly used to model the broad-band continuum emission in AMSPs (Poutanen 2006; Di Salvo & Sanna 2022). The parameters of the thermal Comptonization continuum model are the asymptotic power-law photon index, Γ; the hot, Comptonizing electron temperature, kT_e; and the seed photon temperature, kT_seed. We set the inp_type parameter to zero, indicating that the seed photons are distributed as a blackbody. Since no high-energy cut-off was observed in the XMM-Newton spectrum, we fixed kT_e to 30 keV (see, e.g., Di Salvo et al. 2019). This choice had no significant impact on the obtained results. Using this model, we obtained an unacceptably high reduced χ² of 6 for 146 d.o.f.. The residuals shown in the middle panel of Fig. 9 suggest the presence of a reflection component. We therefore convolved the nthComp component with the rfxconv model from Done & Gierliński (2006) and the relativistic blurring kernel rdblur (Fabian et al. 1989). The parameters of the reflection model are the relative reflection normalization, the iron abundance relative to the solar one, the ionization parameter of the disk’s log ξ, the inner- and outer-disk radii, the disk inclination, and the emissivity index, β, which describes the emissivity of the illuminated disk as a function of the emission radius (∝r^β). We first allowed the iron abundance to vary, but as this did not produce a significant decrease in the fit’s χ², we fixed it at 1. Given the limited spectral coverage and the weakness of the reflection features, we fixed the reflection fraction (rel_refl) at −0.62, consistent with the value found in previous studies (Di Salvo et al. 2019). We also accounted for residual emission features evident in the soft-energy part of the spectrum by introducing three Gaussian components. Some of these features are listed in the official XMM-Newton user guide⁴ as known instrumental lines (the Si-K line at 1.84 keV and the Au-M line at 2.2 keV). We also observed an emission line at 3.5 keV, which we attributed to Ar XVIII fluorescence (see, e.g., Pintore et al. 2016 and Malacaria et al. 2025). In some cases, the spectral resolution of the XMM-Newton/EPIC-pn instrument was insufficient to constrain the width of the modeled lines. Therefore, we fixed the σ parameter to 0 keV. The resulting model gave a reduced χ² of 1.26 with 142 d.o.f. Figure 9 shows the model and the residuals (bottom panel), while Table 2 lists the best-fit parameters.

Fig. 9. XMM-Newton/EPIC-pn spectrum (1.0–10 keV) of high-flux emission from SAX J1808 (top) and best-fit model plotted with a solid line and given in Table 2. The fit model is TBabs*(bbodyrad + gaussian + gaussian + gaussian + nthComp + rdblur*rfxconv*nthComp). The model components are also plotted as dashed lines: the blackbody component, bbodyrad, in red; the Comptonization component, nthcomp, in green; the reflection component in orange; and the Gaussians in purple and yellow. The bottom panel displays the residuals relative to the best-fit model, while the middle panel shows the residuals from a fit performed without the reflection component.

We then modeled the spectrum of the low-flux part of the light curve in the 1–10 keV energy band. We adopted the same spectral model used for the high-flux emission, but without the Gaussian components that described the low-energy emission features, which were not detectable in this spectrum, likely due to lower statistics. Figure 10 shows the model and the residuals, resulting in a reduced χ² of 1.19 with 134 d.o.f. The best-fit parameters are reported in Table 2.

Fig. 10. XMM-Newton/EPIC-pn spectrum (1.0–10 keV) of low part of light curve of SAX J1808 (top) and residuals with respect to best-fit model plotted with a solid line and given in Table 2. The fit model is TBabs*(bbodyrad + nthComp + rdblur*rfxconv*nthComp). The model components are also plotted as dashed lines, with the same color-coding as in Fig. 9.

5. Spectral energy distribution

Figure 11 shows the UV spectrum of SAX J1808, covering the wavelength range of 1650−3100 Å. The only notable feature is a possible weak absorption line near 2600 Å, most likely due to Fe II.

Fig. 11. UV spectrum of HST/STIS observation of SAX J1808 performed on September 10, 2022 and binned at 5 Å in the 1650−3100 Å wavelength range.

We corrected the UV spectrum for interstellar extinction using the empirical relation A_V = N_H/(2.87 ± 0.12)×10²¹ cm⁻² (Foight et al. 2016), where N_H = 2.1 × 10²¹ cm⁻² represents the hydrogen column density along the line of sight to SAX J1808 (Di Salvo et al. 2019), which gives A_V ≃ 0.73 mag. We obtained a corresponding reddening of E(B − V)≃0.27 using the extinction coefficients of Schlafly & Finkbeiner (2011) for R_V = 3.1. We then applied the extinction law of Fitzpatrick (1999) to compute A_λ across the STIS spectral range and corrected each wavelength bin with the corresponding extinction value to derive the dereddened spectrum. We numerically integrated the spectrum with the trapezoidal rule, which approximates the area under the curve by dividing it into small trapezoids, obtaining a UV-flux value of F_UV = (3.7 ± 0.7)×10⁻¹² erg cm⁻² s⁻¹. Considering the measured UV pulse amplitude (see Sect. 3.2), the coherent UV pulses detected during the HST/STIS observation showed a pulsed luminosity of $L_{\mathrm{pulsed}}^{\mathrm{UV}}=0.019L_{\mathrm{UV}}\approx (1.1\pm 0.3)\times {10}^{32}{\mathrm{erg}\mathrm{s}}^{-1}$.

Figure 12 shows the spectral energy distribution (SED) of the total and pulsed emission of SAX J1808 in the UV and X-ray bands during the September 10, 2022 simultaneous observations. To evaluate the X-ray pulsed flux simultaneously with the HST observations, we measured the r.m.s. amplitude in seven energy bands (Fig. 6) and multiplied these values for the unabsorbed X-ray fluxes calculated over the same energy ranges. We then scaled the resulting values by the ratio between the mid-point energy of the interval and the width of the energy interval to convert the pulsed X-ray fluxes into νF_ν units.

Fig. 12. Unabsorbed SED of SAX J1808 from UV to X-rays obtained using simultaneous XMM-Newton/HST (STIS) observations performed during the 2022 outburst reflare phase (September 10, 2022). The HST spectrum is shown in orange from 165 to 310 nm. The pulsed UV flux is plotted with a red triangle. The total 0.5–10 keV X-ray fluxes observed by XMM-Newton are plotted using light blue points. The pulsed X-ray fluxes, plotted as dark blue squares, are calculated over the 0.5–1, 1–2, 2–3, 3–4, 4–5, 5–7, and 7–10 keV energy bands.

6. Discussion

6.1. Phase shifts and hot-spot displacements

We analyzed an uninterrupted XMM-Newton observation of the AMSP SAX J1808 performed during the final reflaring stage of the 2022 outburst. X-ray coherent pulsations were detected throughout the observation, providing strong evidence for channeled accretion onto the magnetic poles of the NS at very low X-ray luminosity. The 0.5–10 keV luminosity measured during the phase of minimum observed X-ray flux (red band in Fig. 2, Table 2) was $L_{\mathrm{X}(\mathrm{low})0.5-10}\simeq 6.{21}_{-0.15}^{+0.20}\times {10}^{34}d{}_{3.5}^2\mathrm{erg}{\mathrm{s}}^{-1}$. This value is on the order of the lowest X-ray luminosity at which pulsations of SAX J1808 were previously observed (around a few ×10³⁴ erg s⁻¹; Hartman et al. 2008, 2009; Patruno et al. 2009a; Bult et al. 2019b; Archibald et al. 2015).

To estimate the location of the inner-disk radius during the observation, we used the relation $R_{\mathrm{m}}\simeq 0.5{\mu }^{4/7}{\dot{M}}^{-2/7}{(2GM)}^{-1/7}$ (Pringle & Rees 1972; Lamb et al. 1973; Campana et al. 2018), where μ = BR_*³ ∼ 10²⁶ G cm³ is the magnetic moment measured from the secular spin-down rate (Hartman et al. 2008), and Ṁ is related to the luminosity as Ṁ = L_X R_∗/GM_∗. We estimated the bolometric 0.1–100.0 keV X-ray luminosity during the low- and high-flux parts of the observations by extrapolating the best-fitting spectral model from the 1.0–10 keV range. We obtained $L_{\mathrm{X}(\mathrm{low})0.1-100}\simeq 0.{91}_{-0.06}^{+0.19}\times {10}^{35}{d}_{3.5}^2{\mathrm{erg}\mathrm{s}}^{-1}$ and $L_{\mathrm{X}(\mathrm{high})0.1-100}\simeq 5.{77}_{-0.12}^{+0.20}\times {10}^{35}{\mathrm{erg}\mathrm{s}}^{-1}$, respectively. The corresponding estimates of the inner-disk radius are $R_{\mathrm{m}(\mathrm{low})}\simeq 36.{55}_{-0.13}^{+0.40}\mathrm{km}$ and $R_{\mathrm{m}(\mathrm{high})}\simeq 21.6_{-0.2}^{+0.4}\text{km}$. Note that these uncertainties are likely underestimated, as they do not include contributions from uncertainties of the magnetic-field strength, mass, and radius. For comparison, the corotation radius of SAX J1808 is $R_{\mathrm{co}}=31m{}_{1.4}^{1/3}\mathrm{km}$, where m_1.4 is the NS mass in units of 1.4 M_⊙. The observation of X-ray pulsations in the low-flux interval confirms that channeled accretion continues even when the mass-accretion rate is low and that the inner disk remains close to the corotation radius, likely in a weak propeller regime (Patruno et al. 2016). Notably, the X-ray pulsations showed a larger amplitude at the lowest observed fluxes (second panel of Fig. 4). This behavior has been observed in several AMSPs, where the pulsed fraction tends to increase as the accretion rate drops (see, e.g., IGR J17379-3747, MAXI J1957+032, MAXI J1816-195, IGR J17498-2921; Bult et al. 2019a, Sanna et al. 2022a, Bult et al. 2022, and Illiano et al. 2024).

A phase-coherent timing analysis revealed a clear phase jump of Δϕ ≃ 0.4 cycles, coincident with a decrease by a factor of three of the X-ray luminosity to the lowest value measured during our observation (L_X(low) ≃ 10³⁵ erg s⁻¹). Timing noise is known to affect the phases observed from SAX J1808 (Hartman et al. 2008, 2009) and other AMSPs (Patruno et al. 2009b). Burderi et al. (2006) observed a phase shift of the opposite sign and roughly half an order of magnitude (Δϕ ≃ −0.2) during the 2002 outburst; this was simultaneous to a rapid drop of the X-ray flux and before the beginning of the reflaring phase. Several models have been proposed to explain such changes. Poutanen et al. (2009) and Ibragimov & Poutanen (2009) interpreted the phase shift and the appearance of a secondary peak at the end of the 2002 outburst in terms of the opening of the view to the antipodal spot due to the recession of the inner disk when the mass-accretion rate decreased. A related interpretation involves changes in the disk-magnetosphere coupling and emission geometry (Kajava et al. 2011). However, in our observation, no secondary peak appeared when the X-ray flux dropped. Instead, the amplitude of the second harmonic became weaker, while the first harmonic became twice as strong compared to the higher flux level (see Figs. 3 and 4). Assuming the drop of the X-ray flux causes an increase in R_in, the lack of a strong secondary peak suggests that the disk was already sufficiently far from the NS to avoid significantly occulting the secondary pole. Indeed, the observed X-ray flux (≲1.5 × 10⁻¹⁰ erg cm⁻² s⁻¹; 0.5–10 keV) was lower than the typical flux at which the transition to the rapid drop occurred in 2002 (≳10⁻⁹ erg cm⁻² s⁻¹). An alternative explanation involves scattering of the hot-spot emission in the accretion funnel (Ahlberg et al. 2024), but this effect is only relevant at higher accretion rates (Ṁ ≳ 10⁻¹⁰ Ṁ yr⁻¹) than those observed in our observation.

Another possible explanation for the observed phase shifts in AMSPs involves spot movement in the azimuthal direction (Lamb et al. 2009; Kulkarni & Romanova 2013). In particular, variations in the hot-spot longitude directly affect the pulse arrival times. Such drifts can result from changes in the mass-accretion rate. As the accretion disk recedes, the inflowing matter couples to different magnetic-field lines and is channeled to slightly displaced regions on the NS surface. The largest phase variation we observed in the XMM-Newton data taken during the reflaring stage of the 2022 outburst corresponds to a shift in the hot-spot longitude of ≃100°. Patruno et al. (2009b) first proposed the existence of a correlation between the observed pulse phase and the X-ray flux, based on the improvement it produced in the description of the phase evolution of a sample of AMSPs. Thus, to account for the phase shifts observed in our analysis, we added a phase-flux correlation term to the constant frequency model, R_flux(t) = b (F_X/F₀)^Λ, as in Bult et al. (2020). This addition yielded a significant improvement in the description of the temporal evolution of the observed pulse phases. The value of Λ = −0.17(9) we obtained agrees with the value suggested by numerical simulations of accretion onto a rapidly rotating NS (Kulkarni & Romanova 2013). These simulations suggest that the azimuthal position of the hot spot relative to the magnetic pole scales proportionally with the magnetospheric radius, δϕ ∝ R_m. Assuming $R_{\mathrm{m}}\propto {\dot{M}}^{\mathrm{\Lambda }}$ and Ṁ to be proportional to the bolometric flux, we obtain δϕ ∝ R_m ∝ F_X^Λ. The value we found is compatible with the index Λ = −1/5, which we found by defining the inner-disk radius as the point where the inflowing matter is forced into corotation with the NS magnetosphere (Spruit & Taam 1993; D’Angelo & Spruit 2010). An even steeper dependence of approximately the usual definition of the Alfvén radius, Λ = −2/7, cannot be excluded, however.

The contour plot in the (b, Λ) plane reveals a correlation between these two parameters, b ⋅ Λ ≃ const (dashed black line in Fig. 5). Such a correlation likely arose because the fit was primarily driven by the step-like variation of Δϕ ≃ 0.4 in response to a variation in the flux Δ F_X ≃ 50 c s⁻¹, with a similar shape. Consequently, the fit was mainly sensitive to the ratio Δϕ/ΔF_X ≃ 8 × 10⁻³ s, which approximates the derivative of the functional form of the phase-flux correlation, dϕ/dF_X = bΛ  ⋅ F₀ ⋅ (F_X/F₀)^Λ − 1 ∼ b ⋅ Λ. Observations of a few phase jumps of variable amplitude, simultaneous to flux variations of different magnitudes, would provide a stronger test of the model and likely remove the degeneracy.

Analysis of NICER observations from the 2019 and 2022 outbursts yielded consistent values of Λ within the uncertainties. Whereas for the 2019 dataset Bult et al. (2020) fixed the Λ index to the theoretically expected value of −1/5, we treated Λ as a free parameter, obtaining Λ = −0.14(5). Our analysis of 2022 data gave a value of Λ = −0.37(14), which differs from that reported by Illiano et al. (2023b) (Λ = −0.81(12)). Such a difference is due to our decision not to consider the phase evolution after the first long gaps in data, as these would have prevented us from reliably approximating the flux with a spline function.

When the accretion disk recedes toward the end of an outburst, the spot latitude is also expected to change as matter is channeled toward the NS surface by magnetic lines that land closer to the magnetic pole. Such variations modify the amplitude of the harmonic components of the pulsed profile (Poutanen & Beloborodov 2006; Lamb et al. 2009). In particular, if the hot spot is close to the stellar spin axis, modest changes in its latitude can produce significant variations. The effect of such changes on the observed pulsed fraction depends on the initial spot latitude (see Fig. 1 from Lamb et al. 2009 and Fig. 6 from Poutanen & Beloborodov 2006). If only a single spot is visible, a decrease in the spot latitude generally reduces the pulsed fraction. On the other hand, if two antipodal spots are visible and the primary spot latitude is ≳45°, its decrease might produce an increase in the fractional amplitude. Simultaneous to the flux decrease and the phase shift discussed above, we observed an increase in the amplitude of the fundamental frequency up to ≃10% (second panel of Fig. 4). Such an increase in the pulsed fraction might thus be interpreted in terms of a change in the spot latitudes caused by a receding inner disk, provided that SAX J1808 is an approximately orthogonal rotator whose pulse profile arises from two spots. In such a geometry, the presence of two nearly antipodal hot spots is also expected to enhance the polarimetric signature. This effect was recently observed by IXPE in the orthogonal rotator GRO J1008−57 (Tsygankov et al. 2023), suggesting that similar measurements of SAX J1808 could provide additional constraints on its geometry.

6.2. UV pulsations

We presented a timing analysis of an HST/STIS UV observation performed simultaneously with the interval of minimum X-ray flux observed by XMM-Newton (Fig. 2). We confirmed the significant detection of UV pulsations, with a pulsed luminosity of $L_{\mathrm{pulsed}}^{\mathrm{UV}}=1.1(3)\times {10}^{32}\mathrm{erg}{\mathrm{s}}^{-1}$ consistent with that observed during the previous accretion event of SAX J1808 in 2019 (Ambrosino et al. 2021).

In both the 2019 and 2022 outbursts, UV pulsations were detected during the reflaring phase (see our Fig. 2 and Fig. 1 of Ambrosino et al. 2021), across an X-ray luminosity range spanning a factor of ∼2 (see Table 3). The values of the pulsed X-ray and UV luminosity we observed simultaneously in the 2022 outburst (see Table 3) give a ratio of $L_{\mathrm{pulsed}}^{\mathrm{X}}/L_{\mathrm{pulsed}}^{\mathrm{UV}}=41\pm 13$. Applying the same analysis to the non-simultaneous 2019 observations yields a similar ratio of 48 ± 18.

Standard mechanisms hardly explain the very bright, pulsed optical and UV flux of SAX J1808. In fact, the observed values imply a brightness temperature of > 1 MeV for a 10 km-wide NS at the distance of SAX J1808, which can be safely excluded given the flux observed in the X-ray band (Ambrosino et al. 2021). On the other hand, the thermal emission observed at soft X-rays is expected to produce a 165–310 nm UV luminosity of

$$\begin{array}{c}\hfill L_{\mathrm{thermal}}=A{\displaystyle {\int }_{{\nu }_1}^{{\nu }_2}}(2\pi k_{\mathrm{B}}{T}_{\mathrm{bb}}{\nu }^2/c^2)\mathrm{d}\nu ,\end{array}$$(3)

where k_B is the Boltzmann constant, kT_bb is the blackbody temperature, A ≈ πR_*²(R_*/R_co) is the hot-spot area of the accreting polar cap (∼100 km² if the whole NS is involved; Frank et al. 2002), and ν₁ and ν₂ represent the lower and upper frequency limits of the UV band. Using kT_bb = 0.362 keV (Table 2), the expected UV luminosity is ∼8 × 10²⁷ erg s⁻¹, which is four orders of magnitude lower than observed. On the other hand, cyclotron emission from hot electrons in the post-shock region of the accretion columns of a 3.5 × 10⁸ G magnetic-field pulsar can account for a UV luminosity up to ≃6 × 10²⁹ erg s⁻¹ (Ambrosino et al. 2021), which is lower by a factor of 150 than the value reported here. Even assuming extreme beaming, it seems highly unlikely that such processes can account for the brightness of the UV pulsations observed from SAX J1808.

The transitional millisecond pulsar PSR J1023+0038 is the only other millisecond pulsar known to show both UV and optical pulsations (Ambrosino et al. 2017; Papitto et al. 2019; Ambrosino et al. 2021; Jaodand et al. 2021; Miraval Zanon et al. 2022). These signals were detected when the source was in an X-ray-faint (L_X ≃ a few × 10³³ erg s⁻¹) accretion-disk regime. The SED of the pulsed flux of PSR J1023+0038 was compatible with a power-law relation connecting the optical and UV to the X-ray band (Papitto et al. 2019; Miraval Zanon et al. 2022). Recently, Baglio et al. (2025) showed that the SED of polarized emission also follows a compatible trend. Synchrotron emission from electrons accelerated where the rotation-powered electromagnetic wind of the pulsar interacts with the mass inflow has been proposed to explain the high pulsed luminosity observed at optical and UV energies (Papitto et al. 2019; Veledina et al. 2019; see also Illiano et al. 2023a; Baglio et al. 2023). Unlike the case of PSR J1023+0038, the simultaneous estimates of the pulsed X-ray and UV fluxes of SAX J1808 cannot be connected by a single power law (see Fig. 12), arguing against an interpretation in terms of a standard synchrotron-emission process. Most importantly, the observed X-ray luminosity of SAX J1808 (L_X ≃ a few × 10³⁵ erg s⁻¹) ensures that the system is accretion powered and that high-density plasma engulfs the magnetosphere. Ambrosino et al. (2021) suggested the possibility of a coexistence or a fast alternation between optical and UV pulsed radiation produced by electrons accelerated by a rotation-powered pulsar and X-ray pulsed emission driven by mass accretion onto the polar caps of the NS. While future theoretical investigation will assess the feasibility of such an unusual scenario, an immediate explanation of the high pulsed luminosity of SAX J1808 that does not rely on strong beaming of the pulsed flux is lacking.

6.3. Evolution in the X-ray spectrum

A study of the evolution of the HR of the source during the XMM-Newton observation revealed a softening of the emission corresponding to the minimum of the source’s X-ray flux. To investigate the origin of this behavior, we performed spectral analysis of both the high- and low-flux intervals of the light curve. The spectra were fit with a model consisting of a blackbody component, a Comptonization component, and a reflection component. Most of the best-fit parameters obtained for the high- and low-flux intervals are consistent at the 3σ confidence level. The normalization of the Comptonization component was the only parameter that decreased significantly in the low-flux interval. The normalization of the bbodyrad component, instead, shows only a slight variation, barely reaching the 3σ confidence level. The ratio between the flux in the Comptonization and blackbody components (see Table 2) changed from $6.4_{-0.4}^{+0.5}$ in the high-flux state to 2.3 ± 0.4 in the low-flux state. Such a behavior could be related to the different physical origins of the two components: the thermal blackbody emission likely originates from a fraction of the NS surface; while the Comptonization emission arises from up-scattering of these same seed photons by electrons in the accretion columns (Poutanen 2006). A decrease in mass accretion rate might then reduce the efficiency of the Comptonization process faster than it affects the thermal blackbody emission, possibly leading to a softer spectrum at lower fluxes.

The parameters measured for the soft and reflection components can provide constraints on the geometry of the emitting regions. The normalization of the bbodyrad component can be used to estimate the radius of the emitting blackbody region, $R_{\mathrm{bb}}=\sqrt{Nor{m}_{\mathrm{bbodyrad}}}·D_{10}$, where D₁₀ = 0.35 is the distance to the source in units of 10 kpc (Galloway & Cumming 2006). We thus obtained $R_{\mathrm{bb}}=3.{26}_{-0.07}^{+0.10}\mathrm{km}$ for the high-flux emission spectrum and R_bb = 2.86 ± 0.14 km for the low-flux part of the light curve. The blackbody component can be interpreted as unscattered emission from a portion of the NS surface. Interestingly, the region appears smaller at lower flux values, possibly suggesting a shrinkage of the emission region. The best-fit parameters of the rdblur convolution model of the high-flux spectrum returned an estimate for the inner-disk radius of R_in < 6.7 R_g (in units of the gravitational Schwarzschild radius R_g = GM/c²), and $i=({58}_{-1}^{+2})°$ for the inclination. The inner-disk radius is equivalent to values of < 14 km for a NS mass of 1.4 M_⊙, which is consistent with the inner disk being truncated before the NS surface and inside the corotation boundary. However, the R_in value inferred from the reflection model appears smaller than expected if the disk is truncated near the corotation radius, as might be suggested by the observed low flux (see Sect. 6.1 for further discussion). Both the inner radius and inclination are consistent with previous estimates based on the relativistic broadening of the reflection features of this source (Papitto et al. 2009; Cackett et al. 2009; Di Salvo et al. 2019). An optical estimate based on quiescent light-curve modeling yields a lower inclination of ${50}_{-5}^{+6}$ degrees (Wang et al. 2013). Such a low inclination would imply an unusually low NS mass to satisfy the companion’s radial-velocity constraints. This would, in turn, exacerbate the tension with the higher inclination required by the X-ray reflection analysis (see Sect. 4.3 of Di Salvo et al. 2019 for further discussion). The uncertainties on the reflection parameters are coarser. This hampers the detection of a change in the inner-disk radius, as suggested by the variation in pulse amplitude and phase.

7. Conclusions

We investigated the temporal and spectral properties of simultaneous X-ray and UV pulsations from the AMSP SAX J1808. We analyzed XMM-Newton and HST observations obtained at low accretion rates during the final reflaring stage of its 2022 outburst. The main results are as follows.

1. We detected significant X-ray pulsations down to a luminosity of $L_{\mathrm{X}0.5-10}\simeq 6.{21}_{-0.15}^{+0.20}\times {10}^{34}{d}_{3.5}^2\mathrm{erg}{\mathrm{s}}^{-1}$ in the 0.5–10 keV band. This value is of the order of the lowest X-ray luminosities at which X-ray pulsations have ever been detected in this source (approximately a few 10³⁴ erg s⁻¹; Bult et al. 2019b).

2. XMM-Newton data revealed significant variations in both pulse phase and amplitude. In particular, a distinct phase jump of approximately 0.4 was observed, coinciding with an amplitude increase up to ∼10% and a decrease in X-ray flux to the lowest luminosity reached during our observation. We interpret these shifts as the drift of emission regions in longitude and latitude on the NS surface as a function of the mass-accretion rate.

3. To take into account the phase shifts, we modeled the phase delays using the flux-adjusted model (Bult et al. 2020), obtaining an index of correlation between the phases and the X-ray flux of Λ = −0.17(9). This value matches the expected scaling of the magnetospheric radius with the accretion rate, which is predicted to be −1/5 from numerical simulations of accretion onto a rotating NS (Kulkarni & Romanova 2013).

4. During the X-ray flux minimum, simultaneous HST observations confirmed significant UV pulsations. The pulsed UV luminosity ($L_{\mathrm{pulsed}}^{\mathrm{UV}}=1.1(3)\times {10}^{32}\mathrm{erg}{\mathrm{s}}^{-1}$) was consistent with that observed during the 2019 outburst. In addition, the pulsed X-ray-to-UV luminosity ratio remained consistent across the two outburst events. The SED of the pulsed emission cannot be explained by a single power-law or thermal model, challenging standard emission scenarios.

5. Spectral analysis of high- and low-flux intervals of the XMM-Newton observation revealed a softening of the emission during the minimum of the X-ray source flux, associated with a decrease in the Comptonization component. The measured inclination ($i=({58}_{-1}^{+2})°$) and the inner-disk radius (< 14 km, for a NS mass of 1.4 M_⊙) remain consistent with previous measurements, supporting a truncated disk within the corotation radius.

Acknowledgments

We thank the referee for useful comments. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA; the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute (program GO/DD-17245), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555; NASA through the NICER mission and the Astrophysics Explorers Program. This work was supported by INAF (Research Grant ‘Uncovering the optical beat of the fastest magnetized neutron stars (FANS)’ and ‘Polarized X-rays from an accreting millisecond pulsar: a pathway to the equation of state of neutron stars (PULSE-X)’, PI: Papitto), the Italian Ministry of University and Research (MUR PRIN 2020) Grant 2020BRP57Z, ‘Gravitational and Electromagnetic-wave Sources in the Universe with current and next-generation detectors (GEMS)’, PI: Astone) and the Fondazione Cariplo/Cassa Depositi e Prestiti (Grant 2023-2560 ‘Taking the optical pulse of the quickest spinning Neutron Stars: a pilot exploratory study (SPES)’, PI: Papitto).

References

All Tables

All Figures

Fig. 1. Light curve of 2022 outburst of SAX J1808 in 0.5–10 keV energy band, using 1 ks bins. NICER observations are shown in black, and XMM-Newton observations are given in red. The blue star marks the epoch of the HST/STIS observation (2022 September 10). We rescaled the XMM-Newton count rate to the NICER count rate using a conversion factor of 1.44, which was calculated with the WebPIMMS tool (https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl), and a power-law model with a photon index of Γ = 2.04 (Table 2).

In the text

Fig. 2. Background-subtracted 0.5–10 keV XMM-Newton/EPIC-pn light curve of 2022 outburst of SAX J1808 measured using 200 s-long bins. A red band indicates the low-flux interval, and a blue band shows the interval of the HST observation discussed in this paper (see Sect. 3.2 for details). Two Type I X-ray bursts, detected at 59831.87 and 59832.82 MJD, are not plotted. Green dash-dotted lines indicate the epochs of their occurrence.

In the text

	Fig. 3. X-ray pulse profiles obtained by folding the high- and low-flux intervals of XMM-Newton observation into 16 phase bins at a spin frequency of νF = 400.975209557 Hz. The solid line represents the best-fit model (constant plus two harmonics), while the dashed red and blue lines indicate the first and second harmonic contributions, respectively. For clarity, two cycles are shown.
In the text

Fig. 4. Top panel: 0.5–10 keV XMM-Newton/EPIC-pn light curve using 1 ks bins, with Type I X-ray bursts removed. The dashed blue line shows the spline interpolation obtained with the UnivariateSpline function from scipy.interpolate. The blue band indicates the interval of the HST observation discussed in this paper (see Sect. 3.2 for details). Second panel: Pulse fractional amplitude for the first harmonic (black dots) and the second harmonic (red dots). Third and fourth panels: Phase residuals relative to a linear model and a model that includes a phase-flux correlation term, respectively. Note that, where not visible, the error bars are smaller than the data points.

In the text

Fig. 5. Contour levels of χ2 obtained by varying both parameters b and Λ in the flux-adjusted phase fit. The contour levels are shown for χmin2 + 2.3 (blue) and χmin2 + 4.61 (red), corresponding to confidence levels of 68% and 90%, respectively, for a fit with two parameters (Lampton et al. 1976). The black dot indicates the values of b and Λ corresponding to the minimum χ2, while the dashed black line represents the contour corresponding to the constant value given by the product b × Λ.

In the text

	Fig. 6. Fractional amplitude of fundamental harmonic of X-ray pulse profiles in seven energy bands (0.5–1, 1–2, 2–3, 3–4, 4–5, 5–7, 7–10 keV) for the low-flux part of the light curve (top panel) and for the high-flux emission (bottom panel).
In the text

	Fig. 7. Pulse profile obtained by folding the whole HST/STIS observation (GO/DD-17245, PI Miraval Zanon) into eight phase bins at a spin frequency of νF = 400.975209557 Hz. The solid line represents the best-fit model, consisting of a single sinusoidal component. For clarity, two cycles are shown.
In the text

Fig. 8. Panel (a): Evolution of HR of SAX J1808 during XMM-Newton observation of reflare phase of 2022 outburst. The HR was estimated from the (2–10 keV)/(0.7–2 keV) energy band. For comparison, the XMM-Newton light curve is shown in blue in the background. Panel (b): Hardness-intensity diagram (HID) of SAX J1808 obtained using XMM-Newton EPIC-pn observation and 200 s bins. The intensity is given in counts per second. The points corresponding to the minimum of the source X-ray flux on the light curve, between 59832.02 MJD and 59832.14 MJD, are marked in red.

In the text

Fig. 9. XMM-Newton/EPIC-pn spectrum (1.0–10 keV) of high-flux emission from SAX J1808 (top) and best-fit model plotted with a solid line and given in Table 2. The fit model is TBabs*(bbodyrad + gaussian + gaussian + gaussian + nthComp + rdblur*rfxconv*nthComp). The model components are also plotted as dashed lines: the blackbody component, bbodyrad, in red; the Comptonization component, nthcomp, in green; the reflection component in orange; and the Gaussians in purple and yellow. The bottom panel displays the residuals relative to the best-fit model, while the middle panel shows the residuals from a fit performed without the reflection component.

In the text

	Fig. 10. XMM-Newton/EPIC-pn spectrum (1.0–10 keV) of low part of light curve of SAX J1808 (top) and residuals with respect to best-fit model plotted with a solid line and given in Table 2. The fit model is TBabs*(bbodyrad + nthComp + rdblur*rfxconv*nthComp). The model components are also plotted as dashed lines, with the same color-coding as in Fig. 9.
In the text

	Fig. 11. UV spectrum of HST/STIS observation of SAX J1808 performed on September 10, 2022 and binned at 5 Å in the 1650−3100 Å wavelength range.
In the text

Fig. 12. Unabsorbed SED of SAX J1808 from UV to X-rays obtained using simultaneous XMM-Newton/HST (STIS) observations performed during the 2022 outburst reflare phase (September 10, 2022). The HST spectrum is shown in orange from 165 to 310 nm. The pulsed UV flux is plotted with a red triangle. The total 0.5–10 keV X-ray fluxes observed by XMM-Newton are plotted using light blue points. The pulsed X-ray fluxes, plotted as dark blue squares, are calculated over the 0.5–1, 1–2, 2–3, 3–4, 4–5, 5–7, and 7–10 keV energy bands.

In the text