© The Authors 2026

1. Introduction

High-mass X-ray binaries (HMXBs) are systems in which a compact object, either a neutron star (NS) or a black hole, accretes matter from a massive companion. Many HMXBs host strongly magnetized NSs (with a magnetic field strength of B ∼ 10¹²–10¹³ G) and, therefore, they often exhibit X-ray pulsations. Among HMXBs, Be/X-ray binaries (BeXRBs) span an exceptionally wide range of accretion rates, providing valuable laboratories for investigating accretion processes, ultrastrong magnetic fields, and quantum electrodynamics (see Mushtukov et al. 2024, for a recent review). In BeXRBs, the accreting matter is supplied through the circumstellar decretion disk around the rapidly rotating Be companion (see Reig 2011, for a review). As a result, these systems exhibit two types of activity: regular type I outbursts occurring near periastron passage, with typical X-ray luminosities of ∼10³⁷ erg s⁻¹ and giant type II outbursts that exceed this level by one or more orders of magnitude.

One such system is the transient source XTE J0111.2−7317, which was discovered by the Rossi X-ray Timing Explorer (RXTE) in November 1998 during a type II outburst, when it reached a luminosity exceeding 10³⁸ erg s⁻¹ (Chakrabarty et al. 1998). Coherent X-ray pulsations with a period of approximately 31 s were detected, confirming its classification as an accreting X-ray pulsar (XRP; Chakrabarty et al. 1998). For convenience, we adopted the short designation SXP31.0, following the naming convention proposed by Coe et al. (2005). SXP31.0 was confirmed to be located in the Small Magellanic Cloud (SMC; Coe et al. 1998; Yokogawa et al. 2000). Simultaneous observations by the Compton Gamma-Ray Observatory (CGRO) detected hard X-ray emission (Wilson & Finger 1998) and follow-up observations with Advanced Satellite for Cosmology and Astrophysics (ASCA) have revealed a pulsating soft X-ray excess. Observations during a quiescent phase in September 2009 showed the source at a luminosity of L_X ≈ 1.8 × 10³⁴ erg s⁻¹ (Christodoulou et al. 2016).

The optical counterpart to SXP31.0 was identified as 2MASS J01110860−7316462 (Coe et al. 2000) and the companion was classified as a B0.5–1Ve star, confirming the system as a BeXRB (Covino et al. 2001; Coe et al. 2003). The orbital period of SXP31.0 was determined to be 90.5 ± 0.1 d (Rajoelimanana et al. 2011; Bird et al. 2012). Early Hα imaging revealed a compact, irregular emission nebula surrounding the star, initially interpreted as a possible supernova remnant (Coe et al. 2000), although spectroscopy would later identify it more plausibly as a locally photoionized H II region (Coe et al. 2003).

To date, observations of SXP31.0 have not yielded a direct measurement of its magnetic field strength, as cyclotron resonant scattering features (CRSFs) have not been detected in its X-ray spectrum (see, e.g., Staubert et al. 2019). Using RXTE data from the 1998 outburst, Kaur et al. (2007) discovered a quasi-periodic oscillation (QPO) at 1.27 Hz and estimated the NS magnetic field strength to be (2–4) × 10¹² G based on the magnetospheric beat frequency model (BFM). The field strength of SXP31.0 thus remains uncertain and can only be inferred through indirect methods.

Recent observations with the Neil Gehrels Swift Observatory (Swift) indicate renewed activity from SXP31.0. The Swift SMC Survey (S-CUBED) detected enhanced X-ray emission from the previously quiescent source on 2025 April 15, with a luminosity of ≈1.8 × 10³⁷ erg s⁻¹ in the 0.3–10 keV band (Gaudin et al. 2025), assuming a distance of 62.44 kpc to the SMC (Graczyk et al. 2020). The source was not detected in the prior survey epoch on March 11, suggesting the onset of a new transient outburst. Subsequently, on 2025 April 26 the Mikhail Pavlinsky ART-XC telescope on board the Spectrum Roentgen Gamma (SRG) observatory detected a bright transient at a position consistent with SXP31.0, with an X-ray flux of (1.2 ± 0.2)×10⁻¹⁰ erg s⁻¹ cm⁻² in the 4–12 keV band (Semena et al. 2025), which corresponds to the luminosity of ∼5 × 10³⁷ erg s⁻¹, consistent with the onset of a Type II outburst.

The Be/X-ray pulsar SXP31.0 has previously exhibited outbursts with luminosities of ∼10³⁸ erg s⁻¹, approaching the Eddington limit, L_Edd for an NS, making it a promising candidate for studying accretion at high rates. Only a handful of super-Eddington pulsars have been identified to date and their properties remain poorly understood. More broadly, XRPs accreting near or above L_Edd provide a valuable opportunity to explore the transition between standard accretion regimes and the radiation-pressure-dominated flows typical of ultraluminous X-ray sources (ULXs), as well as the role of magnetic field interactions in this regime.

In this paper, we present the first comprehensive spectral and timing study of SXP31.0 in the super-Eddington accretion regime, based on observations obtained during its 2025 type II outburst. We triggered dedicated Nuclear Spectroscopic Telescope Array (NuSTAR) and SRG/ART-XC Target of Opportunity (ToO) observations, supported by regular Swift monitoring. Among our key findings is the identification of a previously unreported 0.8 mHz QPO, detected in the early phase of the outburst. The same feature was independently detected by Roy et al. (2025), based on the same NuSTAR ToO dataset. Our preliminary report can be found in Salganik et al. (2025). Here, we present the first detailed spectral and timing analysis and trace the QPO evolution across multiple epochs and luminosity states.

2. Observations and data reduction

Our analysis is based primarily on broadband NuSTAR and SRG/ART-XC observations of SXP31.0, which provide wide energy coverage. These data are supplemented by the Swift data to trace the soft X-ray spectral evolution. The observation log is provided in Table 1.

2.1. NuSTAR observatory

NuSTAR comprises two identical, co-aligned X-ray focal plane modules, FPMA and FPMB (Harrison et al. 2013), which provide sensitive coverage in the 3–79 keV energy band. The NuSTAR observations of SXP31.0 presented here were obtained through a ToO request submitted by our team following the initial identification of the ongoing outburst. Observations were performed on May 2–3, 2025 (ObsID 91101311002; MJD 60797–60798) and May 30–31, 2025 (ObsID 91101314002; MJD 60825–60827), hereafter referred to as NuObs1 and NuObs2, respectively. The source events were extracted using a circular region with a 50″ radius centered on the position reported by Gaudin et al. (2025): RA (J2000) = 01^h11${}^{\mathrm{m}}09\stackrel{\text{s}}{.}51$ (17$\stackrel{°}{.}$78962), Dec (J2000) = $-73°16\prime 44\stackrel{″}{.}9$ ($-73\stackrel{°}{.}27915$), with a 90% confidence error radius of 6.3″. A 100″ radius background region was chosen to optimize the signal-to-noise ratio, particularly at higher energies. The background was extracted from a source-free region in the same detector quadrant (near a detector corner).

Data reduction of the NuSTAR observations was carried out following the standard processing guidelines¹. The data were processed using HEASOFT v6.35.1 and CALDB version 20250415 (including clock correction file 20100101v204). Event files from FPMA and FPMB were barycenter-corrected using the nuproducts tool (via the built-in barycorr option). No significant ghost-ray and stray-light contamination from the nearby bright source SMC X-1 was detected in either FPMA or FPMB for both NuObs1 and NuObs2. Solar activity was elevated during NuObs2, so we filtered outlying events following official recommendations². The spectra and the light curves were extracted with the nuproducts tool as part of the nustardas pipeline. Background-subtracted light curves from FPMA and FPMB were combined using lcmath to improve statistics. No binary motion correction was applied, as the complete set of the source orbital parameters remains unknown.

2.2. SRG/ART-XC

The Mikhail Pavlinsky ART-XC telescope is one of two telescopes installed aboard the SRG observatory, which was launched in 2019 from the Baikonur cosmodrome to the L2 libration point of the Sun-Earth system (Sunyaev et al. 2021). The ART-XC telescope is an imaging instrument consisting of seven modules and operating in photon-counting mode (Pavlinsky et al. 2021) in the 4–30 keV energy band, but the detectors can operate up to ∼120 keV and the highest sensitivity is achieved at 4–12 keV. The time resolution of the instrument is ∼23 μs and the angular resolution is ∼53″. The data were processed using the ARTPRODUCTS v1.0 software with the latest CALDB v20230228. For the spectral and timing analysis, we extracted photons from a circle with a radius of 1$\stackrel{\prime }{.}$8 around the source position and we also applied an energy filtering when needed. Observations performed on 2025 April 30–May 1 (MJD 60795–60796), May 17–18 (MJD 60812–60813), June 6–7 (MJD 60832–60833), and June 19–21 (MJD 60845–60847) are referred to below as ArtObs1, ArtObs2, ArtObs3, and ArtObs4, respectively. The background was estimated by collecting photons from a circular region with a radius of 5$\stackrel{\prime }{.}$4 on the detectors that did not overlap with the region in which the source photons were extracted. SXP31.0 was the only detected source within the ART-XC field of view during the observation.

2.3. Swift observatory

To track the evolution of the outburst (see Fig. 1), we used data from the monitoring campaign with the XRT telescope on board the Swift observatory (Burrows et al. 2005; Gehrels et al. 2004), comprising 35 observations (see Table 1). These observations span MJD 60781.5–60930.5 and were carried out in photon counting (PC) or windowed timing (WT) modes, chosen according to the source count rate. The source and background spectra for each individual observation were extracted using the UK Swift Science Data Centre’s XRT data-analysis software³, which also automatically selected the extraction regions following the method described in Evans et al. (2009).

Fig. 1. Evolution of SXP31.0 during the 2025 outburst based on Swift/XRT, NuSTAR, and SRG/ART-XC observations. The light curve in the 4–12 keV energy range is shown at the top panel. The gray points show the CGRO/BATSE 20–50 keV pulsed luminosity from the 1998 outburst, overplotted for comparison. Along the time axis it was placed arbitrarily, and the luminosity was scaled by a factor of 0.45. The spin period (see Table 2) is shown at the bottom panel.

2.4. Spectral data approximation

To enable a joint spectral fitting, the broadband spectra were grouped into five distinct epochs based primarily on similar flux levels. Observations within each group were taken at comparable luminosities and, coincidentally, they also exhibit consistent pulse profile morphology (see Sects. 3.1.2 and 3.2). Within each epoch, data from NuSTAR, SRG/ART-XC, and Swift/XRT were combined. All spectra were binned with a minimum of one count per energy bin and fitted using W-statistics (Wachter et al. 1979).

Unless stated otherwise, all uncertainties correspond to the 1σ confidence level, and all fluxes reported in this work are unabsorbed. Luminosities were computed assuming a distance of 62.44 kpc to the SMC (Graczyk et al. 2020).

To construct the long-term light curve, individual Swift/XRT spectra were fitted with an absorbed power law plus blackbody model, tbabs × (po + bbodyrad), with the blackbody temperature fixed at T = 0.22 keV (see Sect. 3.3). Throughout this analysis, the hydrogen column density was fixed at N_H = 0.18 × 10²² cm⁻², following Yokogawa et al. (2000). The resulting long-term light curve is shown in Fig. 1. For comparison, the CGRO/BATSE 20–50 keV pulsed flux from the 1998 outburst⁴ (Yokogawa et al. 2000) is also overplotted. For the conversion from Crab units to the energy flux in 20–50 keV range, we adopted 1 Crab = 9.22 × 10⁻⁹ erg s⁻¹ cm⁻² (Frontera et al. 2007). The BATSE light curve of the 1998 outburst differs in shape from the 2025 event: while both reached comparable peak luminosities of ∼10³⁸ erg s⁻¹, the 1998 outburst exhibited a pronounced two-humped shape and its decay was noticeably steeper. However, these comparisons should be taken with caution since the BATSE measurements represent pulsed flux.

3. Results

We performed a detailed timing analysis, including energy-resolved pulse profiles and pulsed fraction measurements, and we report the discovery and modeling of a transient 0.8 mHz QPO. We also present results of the broadband spectral fitting. To investigate the evolution throughout the outburst, the observations were divided into two main intervals: a QPO-active phase (epoch 1; see Table 2), when the 0.8 mHz QPO was detected, and a subsequent QPO-free phase (epochs 2–5).

3.1. Timing analysis during QPO-active interval

3.1.1. Detection of the 0.8 mHz QPO

Upon closer inspection, the light curve of SXP31.0 reveals quasi-periodic variability with a characteristic timescale of ∼10³ s (bottom panel of Fig. 2). To investigate this behavior, we computed a dynamic power density spectrum (PDS) from the NuObs1 barycenter-corrected 3–79 keV light curve using the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982), as implemented in the astropy Python library. Each time bin in the dynamic PDS corresponds to a power spectrum computed over the following 8192-s interval, with a step size of 512 s. Each individual power spectrum was logarithmically rebinned in frequency to enhance visibility (top panel of Fig. 2). The quasi-periodic variability is most prominent during the first 32.8 ks of the observation (MJD 60797.8239–60798.2041). The hatched region on the right-hand side of Fig. 2 marks the interval where the sliding window cannot be applied without extending beyond the observational data span.

Fig. 2. Smoothed dynamic power spectrum (top) and light curve (bottom) of NuObs1 in the 3–79 keV energy band. The dynamic power spectrum is computed using a sliding window of 8192 s with a step of 512 s. The hatched region indicates the part of the time range where the sliding window cannot be applied without extending beyond the data boundaries. The light curve is binned with a 10 s time resolution.

To characterize the QPO, we extracted this segment of the light curve and computed an averaged PDS using the powspec tool from the XRONOS package. The PDS was calculated using root-mean-square (rms) normalization to express variability in terms of fractional rms amplitude. The white noise level was subtracted and the spectrum was also logarithmically rebinned geometrically with a factor of 1.1. The same procedure is used for all other power spectra presented in this work unless stated otherwise.

The resulting PDS shown in Fig. 3 reveals a broad feature near ∼10⁻³ Hz, consistent with a mHz QPO. This QPO is accompanied by an additional excess at higher frequencies, hereafter referred to as the “shoulder.” A similar QPO-shoulder structure, along with low-frequency red noise, was previously observed in the super-Eddington pulsar M51 ULX-7 (Imbrogno et al. 2024). Following the same modeling approach, we fit the PDS using components for the QPO, the shoulder, and the pulsar spin frequency. We excluded the red noise component because the available data do not provide sufficient coverage at low frequencies (below 10⁻⁴ Hz) to constrain it reliably.

Fig. 3. Rms-normalized PDS of SXP31.0 during the QPO-active interval for NuObs1 (3–79 keV) computed from the first 32.8 ks of the observation where the QPO is most prominent. The best-fit model is shown in red, with the green dashed line for the QPO component and the blue dashed line for the shoulder component. Panel (b) shows residuals to the model with only the QPO and spin frequency, while panel (c) shows residuals after adding the shoulder component.

We converted the observed PDS into an XSPEC-compatible format using the flx2xsp utility and fitted it with a sum of three lorentz models in XSPEC, using a Whittle likelihood for the parameter estimation. This approach is more suitable than least squares when dealing with exponentially distributed periodogram powers (e.g., Whittle 1953, 1957). All model parameters were left free during fitting, except for the spin frequency, which was fixed at ν_spin = 1/30.44970 Hz. The best-fit model yielded a QPO centered at ν_QPO = 0.8 ± 0.1 mHz, with a quality factor of $Q=3.2_{-1.8}^{+10.6}$ and a fractional rms amplitude of 14 ± 5%. The shoulder component was located at ν_shoulder ≈ 1.5 mHz, with a poorly constrained quality factor. The best-fit values and their uncertainties are listed in Table 3.

To assess the significance of the detected QPO, we performed Monte Carlo simulations using the simftest⁵ procedure from the XSPEC package, providing the difference in Whittle statistic between the models with and without the QPO component. We ran 10⁴ simulations including the narrow Lorentzian QPO component. The resulting detection significance corresponds to at least 3σ. Thus, the QPO is firmly detected in the course of the NuObs1 observation, when the source was accreting at a bolometric luminosity of L_bol ≈ 2.5 × 10³⁸ erg s⁻¹ (see Table 4).

To verify the presence and stability of the QPO, we independently analyzed the 4–25 keV ArtObs1 data, obtained just a few days earlier (see Fig. 4). The ART-XC PDS also exhibits a QPO at approximately 0.9 mHz (see Table 3), with parameters that are consistent with those observed in NuObs1. For ArtObs1, we included the red noise in the model as a simple power law and the shoulder component was not required for the fit.

Fig. 4. Same as Fig. 3 but for ArtObs1 (4–25 keV) observation. Here the orange dashed line shows the red-noise component. Panel (b) shows the residuals for the model with the red noise only, while panel (c) shows residuals after adding the QPO component.

3.1.2. Pulse profiles

A standard epoch-folding technique, implemented in the efsearch tool of the FTOOLS package, was used to estimate the spin period values for the observations presented in this paper (see Table 2 and the bottom panel of Fig. 1). The uncertainty was estimated from the distribution of measured periods obtained from a set of 10³ simulated light curves, following the procedure from Boldin et al. (2013).

Pulse profiles are a valuable diagnostic tool for investigating the emission characteristics of an NS and the properties of its surrounding environment. For SXP31.0, we generated pulse profiles by folding the barycenter-corrected light curves with the best-fit spin periods (see Table 2) using the efold task from the XRONOS package. The energy-resolved profiles for all epochs are presented in Fig. 5.

The pulse profile during epoch 1 exhibits a three-peaked structure at low energies, which gradually evolves with increasing energy (see Figs. 5a–c). As the energy increases, the secondary peaks weaken: the profile first transitions to a double-peaked shape and then to a single-peaked, nearly sinusoidal form above ∼26 keV. Notably, at these highest energies, a pronounced shoulder appears on the main peak around phases ∼0.2–0.5 (Figs. 5c).

Fig. 5. Energy-resolved pulse profiles of SXP31.0. Gray dashed vertical lines mark the boundaries of the pulse peaks. The phase corresponding to the maximum of the profile is set to zero. Two full periods are shown per profile, each vertically offset for clarity. The flux in each energy band is normalized to its mean value. Energy ranges (in keV) are labeled next to each profile.

To explore these variations in more detail, we analyzed energy-resolved folded light curves to derive the energy dependence of the pulsed fraction (PF), which is commonly defined as [max(rate) – min(rate)]/[max(rate) + min(rate)]. The PF, calculated using pulse profiles that were divided into 15 phase bins, displays a behavior that is unusual for XRPs (as shown in Fig. 6). While most bright XRPs exhibit a monotonic increase in the PF with energy (Lutovinov & Tsygankov 2009), SXP31.0 behaves differently: it remains at the level of ≈10% up to 20 keV with a local minimum near the 6.4 keV iron line. At higher energies, PF begins to rise, reaching values of approximately 35% in the 36–79 keV energy band. This increase in the PF coincides with the previously noted phase broadening of the main peak and the disappearance of the secondary peaks.

Fig. 6. Dependence of the PF on energy. The dashed red line marks the energy above which secondary peaks disappear from the pulse profile in NuObs1, coinciding with the onset of a growth in the PF (see Sect. 3.1.2). L38 is the bolometric luminosity in units of 1038 erg s−1.

3.2. Timing analysis during QPO-free interval

On MJD 60813, epoch 2 observations were performed at a higher luminosity of L_bol = 3.6 × 10³⁸ erg s⁻¹ (see Table 4), which is approximately 1.4 times higher than in epoch 1. One of the most notable results is the complete disappearance of the QPO from the PDS (see Fig. 7a), despite improved statistics owing to the longer exposure and higher flux. Thus, the QPO vanished within 14 days of its last detection. The PDS can be described as a red noise (a simple power law) and the QPO component was not required. However, this transition is accompanied by a sharp increase in PF in the 4–25 keV band, rising from an average value of ≈10% in epoch 1 to ≈26% in epoch 2 (see Fig. 6).

Fig. 7. Rms-normalized PDS of SXP31.0 obtained after the initial detection. Upper subpanels show the data (black crosses) and the best-fit model (red line). The lower subpanels show the residuals. Panels (a), (b), (c), and (d) correspond to ArtObs2, NuObs2, ArtObs3, and ArtObs4 datasets, respectively. See text for details.

An energy-resolved analysis of the pulse profile reveals a similar triple-peaked shape at all energies up to 25 keV (Fig. 5d). However, this structure appears significantly distorted compared to the previously observed pattern at epoch 1, despite only a modest increase in the flux. The relative phase spans of the peaks have changed substantially, indicating a change in the pulse morphology. The associated increase in PF, which remains at ∼26% in all energy ranges up to 25 keV (Fig. 6), along with the altered pulse shape, suggests a change in the beam configuration.

Subsequently, observations of epochs 3 and 4 were performed at a bolometric luminosity of L_bol= (2.7–2.8)×10³⁸ erg s⁻¹ (see Sect. 3.3), which is very similar to that of epoch 1. Despite that, no evidence of the previously detected QPO was found in the PDS (Figs. 7b,c) and the timing properties differed markedly from those observed in epoch 1 (see Fig. 5e,f). The PF remains at approximately 25% up to 20 keV and then gradually increases to about 40% at higher energies (Fig. 6). The 3–79 keV pulse profile of SXP31.0 exhibits a double-peaked structure (Figs. 5e–g). An exception is seen in the 36–79 keV range, where the profile becomes distinctly single-peaked.

Finally, epoch 5 is of particular interest, as it is the only observation obtained at a luminosity lower than that of epoch 1, where the QPO was detected (with epoch 5 reaching L_bol = 2.1 × 10³⁸ erg s⁻¹), yet it still shows no evidence of the QPO (Fig. 7d). The energy-resolved pulse profile retains a double-peaked structure at lower energies (see Fig. 5h); the second peak vanishes above 16 keV. The PF remains moderately higher than in epoch 1, though lower than in the other QPO-free states (epochs 2 to 4), reaching 17% in the 4–25 keV range.

3.3. Spectral analysis

We performed the pulse phase-averaged spectral analysis of SXP31.0 using data from epochs 1–5 in order to trace the evolution of the spectral parameters in different states. We excluded the 3–4 keV range from the NuSTAR data analysis due to cross-calibration difficulties with Swift, and energies below 0.8 and 0.5 keV were ignored for Swift in WT and PC modes, respectively.

To approximate the source spectrum in epoch 1, we adopted the canonical accreting XRP model (e.g., Filippova et al. 2005; Coburn et al. 2002), consisting of a power-law continuum with an exponential cutoff at high energies. However, a simple cutoffpl model fails to reproduce the continuum adequately, leaving structured residuals (see Fig. 8b). To improve the fit, we added a blackbody component with a temperature of T_bb ≈ 1.3 keV, motivated by a similar approach used in modeling other super-Eddington pulsars, such as Swift J0243.6+6124 (Tao et al. 2019; Bykov et al. 2022) and SMC X-3 (Pottschmidt et al. 2016; Tsygankov et al. 2017), where comparable blackbody temperatures were found. In our case, this significantly improved the fit (W-stat./d.o.f. = 1.08, see Fig. 8c). Furthermore, after taking into account the soft excess with T_excess = 0.22 keV and R_excess = 160 km, previously reported in Yokogawa et al. (2000), we were able to reduce W-stat./d.o.f. ratio to the value of 1.01 (Fig. 8d). The inclusion of the hot bbodyrad component is also required when adopting an alternative continuum model such as po×highecut (Coburn et al. 2002), reinforcing its necessity. We were unable to detect the presence of any cyclotron absorption features. The best-fit model parameters are presented in Table 4.

Fig. 8. Unfolded energy spectra of SXP31.0. The spectra were spaced along the y-axis for visual convenience with the factors shown. The black, green, magenta, blue, and red points correspond to the epoch 1–5 data and their respective model approximations. The dashed lines represent the individual model components (cutoffpl, bbodyrad, gauss, and a second bbodyrad) for the epoch 5 dataset. Panels (b), (c), and (d) show the residuals of the epoch 1 dataset when fitted with the cutoffpl model alone, with an additional bbodyrad component, and with two bbodyrad components, respectively.

The spectrum also shows an iron fluorescence line at 6.4 keV, previously reported by Yokogawa et al. (2000), which is modeled with a narrow Gaussian component gauss. The line width σ was fixed at 0.1 keV for the SRG/ART-XC data, while it was left as a free parameter for the NuSTAR observations. A cross-calibration constant is included to account for possible normalization differences between the FPMA, FPMB, ART-XC, and Swift/XRT.

To investigate the changes in the spectral properties, we analyzed the epoch 2 spectra using the same model as in epoch 1, which remained statistically preferred. For epochs 2–5, the soft excess temperature was fixed at T_excess = 0.22 keV to reduce parameter degeneracy, adopting the best-fit value from epoch 1 (see Table 4).

Although ArtObs2 is separated by approximately five days from the closest available Swift/XRT observation (ObsID 00019718011; see Table 1), both datasets exhibit consistent flux levels in the 4–10 keV band, with values of 1.4 × 10⁻¹⁰ erg s⁻¹ cm⁻². This agreement justifies their joint spectral fitting, which enables tighter constraints on spectral parameters. The spectrum was approximated using the same model as in the earlier analysis, and the fit remains statistically acceptable, with a W-stat./d.o.f. ratio of approximately 1.07. We also analyzed the spectra from epochs 3, 4, and 5 using the same model, which provided statistically acceptable fits with W-stat./d.o.f. values of 0.98, 1.14, and 1.11, respectively (Table 4). None of the spectra from any epoch show evidence of CRSFs in the 5–50 keV energy range.

We also analyzed pulse-phase-resolved spectra for NuObs1 and NuObs2, obtained by dividing the data into five equally spaced pulse phase bins. Each spectrum was fitted using the same model as the phase-averaged spectrum except soft excess: const × tbabs × (cutoffpl + bbodyrad + gauss). Swift/XRT data were not included. The energy of the iron line was frozen to the value E_iron = 6.4 keV and the width values were frozen to σ_iron = 0.35 and 0.5 keV for NuObs1 and NuObs2, respectively, based on the best-fit values from the phase-averaged spectral analysis (see Table 4). The results can be found in Fig. 9. We were unable to detect any cyclotron absorption features in any of the spectra. Due to the low counting statistics and complexity of the spectral shape, the significance of the variations in the spectral parameters was hard to assess.

Fig. 9. Dependence of the hardness ratios and the spectral parameters on the pulse phase for the NuObs1 (left) and NuObs2 (right). The averaged pulse profile in a wide energy range 3–79 keV is superimposed in gray for visual comparison.

4. Discussion and conclusions

The NuSTAR, Swift/XRT, and SRG/ART-XC observations of SXP31.0 during its 2025 type II outburst provide the first comprehensive broadband spectral and timing characterization of this XRP. The data span the 0.5–79 keV energy range and enable a detailed investigation of both pulse-phase-averaged and pulse-phase-resolved properties. The source exceeded the Eddington limit of 1.8 × 10³⁸ erg s⁻¹ for a canonical 1.4 M_⊙ NS during the outburst. This firmly places SXP31.0 in the super-Eddington accretion regime, among the most luminous outbursts ever observed in BeXRB systems.

Throughout all broadband observations, the phase-averaged X-ray spectrum was characterized by a hard power-law continuum, modified by a high-energy exponential cutoff. The inclusion of two blackbody components with T ≈ 1.3 keV, R ≈ 7 km, and a cooler excess component with T_excess ≈ 0.2 keV, R ≈ 140 − 170 km was required to accurately reproduce the observed spectrum. A narrow Fe Kα emission line was detected at 6.4 keV, with an equivalent width not exceeding ≈110 eV. No CRSF, a direct probe of the NS’s magnetic field (Staubert et al. 2019), was detected in either phase-averaged or phase-resolved spectra of SXP31.0 across the 5–50 keV range in any of the observations presented in this study, suggesting magnetic field strength B < 5 × 10¹¹ G or > 5 × 10¹² G. A shallow or complex-shaped feature below the detection threshold cannot be entirely ruled out.

We discovered the presence of a 0.8-mHz QPO in SXP31.0 at a luminosity of L_bol ≈ 2.5 × 10³⁸ erg s⁻¹ in the epoch 1 observations. It establishes SXP31.0 as the fourth known super-Eddington XRP to exhibit mHz low-frequency quasi-periodic variability (for a review, see Čemeljić et al. 2025; Veresvarska et al. 2025), and the first such detection in a super-Eddington luminosity state outside the class of confirmed ultraluminous X-ray pulsars (ULXPs). This newly detected QPO does not match the 1.27 Hz feature previously reported by Kaur et al. (2007), which was observed during an earlier outburst at a similar luminosity level of ∼10³⁸ erg s⁻¹.

In NuObs1, the energy-resolved pulse profile of SXP31.0 evolved from a complex three-peaked structure at soft X-rays to a single, nearly sinusoidal shape above 26 keV (Fig. 5b,c). At low energies (≲10 keV), the PF in SXP31.0 remains relatively low, at around ∼10%, consistent with values observed in other super-Eddington XRPs exhibiting mHz QPOs. For example, M51 ULX-7 (P_spin ≈ 2.8 s, ν_QPO ≈ 0.5–0.6 mHz; Imbrogno et al. 2024) shows an upper limit of 7%, NGC 7793 P13 (P_spin ≈ 0.4 s, ν_QPO ≈ 10 mHz; Imbrogno 2024; Israel et al. 2025) exhibits PFs in the range 5–10%, and in case of M82 X-2 (P_spin ≈ 1.3 s, ν_QPO ≈ 2.8–4.0 mHz; Feng et al. 2010), no detectable pulsations were reported during the only known QPO episode. All of these three systems had their QPOs observed at luminosities L_X ≳ 10³⁹ erg s⁻¹.

However, unlike these sources, SXP31.0 was also observed at higher energies, where it exhibits a marked increase in the PF, rising to ∼35% above 30 keV (see Sect. 3.1.2). Since the other systems lack observations above 10 keV, any potential increase in PF at higher energies may have gone unnoticed due to a limited energy coverage. This suggests that hard X-ray observations could help identify NSs in ULXs showing mHz QPOs, even in the absence of low-energy pulsations.

In the subsequent epoch 2 observations at the luminosity of L_bol ≈ 3.6 × 10³⁸ erg s⁻¹ (see Table 4), the QPO was no longer detected. At the same time, the PF in the 4–25 keV energy band increased sharply from approximately 10% to 26%. The pulse-profile shape also exhibited a three-peaked structure, although the overall morphology appeared distorted.

In epochs 3 and 4, observations at the luminosities L_bol≈(2.7–2.8)×10³⁸ erg s⁻¹ indicate that the QPO was once again not detected and the PF showed similar to the epoch 1 behavior, rising to 40% above 25 keV. Energy-resolved pulse profiles evolved into a two-peaked shape, transitioning to a sinusoidal, single-peaked profile above 36 keV.

In epoch 5, at a lower luminosity of L_bol = 2.1 × 10³⁸erg s⁻¹, no QPO was detected. The pulse profile was double-peaked at low energies, with the second peak vanishing above 16 keV, and the PF reached 17% in the 4–25 keV range, somewhat below other QPO-free epochs.

These findings indicate that the QPO exhibits a distinctly transient nature, appearing only under specific physical conditions, since it is absent at both higher and lower luminosities after the detection. Moreover, QPO presence is associated with a marked decrease in the PF.

The drastic frequency shift of more than three orders of magnitude between the 1.27 Hz QPO reported by Kaur et al. (2007) and the 0.8 mHz QPO detected in this work cannot be explained by standard QPO models such as the BFM (Alpar & Shaham 1985) or the Keplerian frequency model (van der Klis et al. 1987), which predict ν_QPO ∝ Ṁ^3/7. Explaining such a shift would require an unobserved change in the accretion rate by nearly seven orders of magnitude. This strongly suggests a different physical origin for the newly detected mHz QPO.

One explanation comes from the disk precession model proposed for ULXPs by Čemeljić et al. (2025), in which global, rigid-body-like precession of a tilted inner accretion disk is driven by magnetic torques from the NS’s inclined dipole. A key feature of this model is the inverse scaling of the QPO frequency with a spin period, ν_QPO ∝ 1/P_spin, which is largely independent of magnetic field strength or accretion rate, and primarily governed by the disk geometry. Specifically, the model requires the presence of a geometrically thick inner accretion disk capable of supporting global precession.

For SXP31.0, with a spin period of 30.45 s, the model predicts a QPO frequency of about 0.1 mHz, placing the observed 0.8 mHz signal within a factor of eight of the expected value. Considering the model’s sensitivity to uncertain geometric parameters, this discrepancy is not necessarily problematic. The proportionality constant, K, in the ν_QPO = K/P_spin relation depends sensitively on the disk geometry and may deviate from the established K ≈ 3 mHz s if the accretion flow in SXP31.0 is not fully radiation-pressure dominated or geometrically thick, as assumed for confirmed ULXPs. Indeed, SXP31.0 resides near the Eddington limit and may not sustain a thick, radiation-pressure-dominated disk typical of confirmed ULXPs (Mushtukov et al. 2019). Instead, it may represent a transitional regime, where the disk is only moderately inflated and geometrically intermediate. In such a configuration, the conditions required for coherent global precession may be only marginally satisfied, potentially shifting the QPO frequency away from the expected scaling.

Alternatively, the observed 0.8-mHz QPO in SXP31.0 may be interpreted within the magnetically driven precession (MDP) framework originally developed by Shirakawa & Lai (2002a,b), a more general version of the global precession scenario later applied by Čemeljić et al. (2025). While the latter relies on a simplified scaling for a radiation-pressure-supported thick disk in spin equilibrium, the MDP model derives the precession frequency from first principles, balancing magnetic torques and viscous stresses in a warped accretion flow. Although the model in principle allows one to constrain the dipolar magnetic field, B, we avoid quoting specific values due to significant uncertainties in key parameters such as the viscosity parameter, α, magnetic obliquity. θ, and the disk thickness factor, D(r). Importantly, the MDP framework of Čemeljić et al. (2025) was calibrated specifically for the ULXP systems, which display mHz QPO, which is an order of magnitude brighter than SXP31.0.

The MDP-based interpretation supports the scenario in which mHz QPOs originate in regions of the accretion flow that also influence the formation or visibility of pulsations. It also underscores the need for caution when using QPO frequencies in ULXs as diagnostics of the compact object mass, since the physical origin of the QPO may not be directly tied to Keplerian timescales. Finally, the apparent anticorrelation between QPOs and PF may partly explain the difficulty in detecting pulsations in many ULXs and suggests that the fraction of NSs among the ULX population could be higher than currently inferred.

Moreover, the current data provide the first observational constraint on how quickly the disappearance of the mHz QPO and the accompanying rise in PF can occur. In SXP31.0, this transition took place within a ≲14 day interval between May 3 and May 18, setting an upper limit on the temporal evolution of this phenomenon in the super-Eddington accretion regime.

Acknowledgments

AS acknowledges support from the EDUFI Fellowship and Jenny and Antti Wihuri Foundation (grant no. 00240331). SST acknowledges financial support from the Centre for Astrophysics and Relativity at Dublin City University during his visit. AAM acknowledges UKRI Stephen Hawking fellowship. SVM, IYL, AAL and AYT acknowledge support from the Ministry of Science and Higher Education of RF grant 075-15-2024-647. We thank the anonymous referee for their careful reading of the paper and thoughtful suggestions, which helped improve the clarity of the results. AS also thanks Vadim Kravtsov for helpful discussions during the preparation of this work. This work is partially based on observations with the Mikhail Pavlinsky ART-XC telescope, hard X-ray instrument on board the SRG observatory. The SRG observatory was created by Roskosmos in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI) in the framework of the Russian Federal Space Program, with the participation of Germany. The ART-XC team thanks the Roscosmos State Corporation, the Russian Academy of Sciences, and Rosatom State Corporation for supporting the ART-XC telescope, as well as the JSC Lavochkin Association and partners for manufacturing and running the Navigator spacecraft and platform. This work made also use of data supplied by the UK Swift Science Data Centre at the University of Leicester and data obtained with NuSTAR mission, a project led by Caltech, funded by NASA, and managed by JPL. This research also has made use of the NuSTAR Data Analysis Software (NUSTARDAS) jointly developed by the ASI Science Data Centre (ASDC, Italy) and Caltech. We are grateful to the Swift and NuSTAR teams for approving and rapid scheduling of the monitoring campaign. This work made use of Astropy: (http://www.astropy.org) a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022). This research has made use of data and software provided by the High Energy Astrophysics Science Archive Research Centre (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory.

References

All Tables

All Figures

Fig. 1. Evolution of SXP31.0 during the 2025 outburst based on Swift/XRT, NuSTAR, and SRG/ART-XC observations. The light curve in the 4–12 keV energy range is shown at the top panel. The gray points show the CGRO/BATSE 20–50 keV pulsed luminosity from the 1998 outburst, overplotted for comparison. Along the time axis it was placed arbitrarily, and the luminosity was scaled by a factor of 0.45. The spin period (see Table 2) is shown at the bottom panel.

In the text

Fig. 2. Smoothed dynamic power spectrum (top) and light curve (bottom) of NuObs1 in the 3–79 keV energy band. The dynamic power spectrum is computed using a sliding window of 8192 s with a step of 512 s. The hatched region indicates the part of the time range where the sliding window cannot be applied without extending beyond the data boundaries. The light curve is binned with a 10 s time resolution.

In the text

Fig. 3. Rms-normalized PDS of SXP31.0 during the QPO-active interval for NuObs1 (3–79 keV) computed from the first 32.8 ks of the observation where the QPO is most prominent. The best-fit model is shown in red, with the green dashed line for the QPO component and the blue dashed line for the shoulder component. Panel (b) shows residuals to the model with only the QPO and spin frequency, while panel (c) shows residuals after adding the shoulder component.

In the text

	Fig. 4. Same as Fig. 3 but for ArtObs1 (4–25 keV) observation. Here the orange dashed line shows the red-noise component. Panel (b) shows the residuals for the model with the red noise only, while panel (c) shows residuals after adding the QPO component.
In the text

	Fig. 5. Energy-resolved pulse profiles of SXP31.0. Gray dashed vertical lines mark the boundaries of the pulse peaks. The phase corresponding to the maximum of the profile is set to zero. Two full periods are shown per profile, each vertically offset for clarity. The flux in each energy band is normalized to its mean value. Energy ranges (in keV) are labeled next to each profile.
In the text

	Fig. 6. Dependence of the PF on energy. The dashed red line marks the energy above which secondary peaks disappear from the pulse profile in NuObs1, coinciding with the onset of a growth in the PF (see Sect. 3.1.2). L38 is the bolometric luminosity in units of 1038 erg s−1.
In the text

	Fig. 7. Rms-normalized PDS of SXP31.0 obtained after the initial detection. Upper subpanels show the data (black crosses) and the best-fit model (red line). The lower subpanels show the residuals. Panels (a), (b), (c), and (d) correspond to ArtObs2, NuObs2, ArtObs3, and ArtObs4 datasets, respectively. See text for details.
In the text

Fig. 8. Unfolded energy spectra of SXP31.0. The spectra were spaced along the y-axis for visual convenience with the factors shown. The black, green, magenta, blue, and red points correspond to the epoch 1–5 data and their respective model approximations. The dashed lines represent the individual model components (cutoffpl, bbodyrad, gauss, and a second bbodyrad) for the epoch 5 dataset. Panels (b), (c), and (d) show the residuals of the epoch 1 dataset when fitted with the cutoffpl model alone, with an additional bbodyrad component, and with two bbodyrad components, respectively.

In the text

	Fig. 9. Dependence of the hardness ratios and the spectral parameters on the pulse phase for the NuObs1 (left) and NuObs2 (right). The averaged pulse profile in a wide energy range 3–79 keV is superimposed in gray for visual comparison.
In the text