Sidebands to mHz QPOs in 4U 1626−67 in the second spin-down state

ABSTRACT

We report results from an AstroSat target-of-opportunity (ToO) observation of 4U 1626−67, performed on 2023 May 18, soon after the discovery of torque reversal to spin-down in the source. The X-ray emission exhibited significant dependence on both energy and torque state. This work highlights the comparison of timing features of 4U 1626−67 with a previous AstroSat observation from 2018, when the neutron star was in the spin-up state. The power density spectrum (PDS) of the 2023 observation comprised a sharp peak corresponding to ${\nu }_{\mathrm{N}\mathrm{S}}\sim$130 mHz X-ray pulsations along with a prominent quasi-periodic oscillation (QPO) feature at ${\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}\sim$46 mHz with $\sim$20 per cent rms amplitude, which was positively correlated with energy. We also report the detection of sidebands to QPO occurring at a beat frequency (⁠${\nu }_{\mathrm{N}\mathrm{S}}-{\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}$⁠) of $\sim$83 mHz with $\sim$8 per cent rms amplitude, having $>3\sigma$ detection significance. Additionally, we utilized Nuclear Spectroscopic Telescope ARray (NuSTAR) observations from the same torque state (2023 May–July) to analogize the presence and energy dependence of sidebands. The source retains timing properties in this spin-down torque state, similar to those seen in the previous spin-down phase. In sharp contrast, PDS from the 2018 observation was dominated by red noise, an absence of QPOs and a broadening in the wings of the pulse frequency peak, indicating a coupling between periodic and low-frequency aperiodic variability. Furthermore, we detected the known cyclotron resonance scattering feature (CRSF) at 37 keV in the Large Area X-ray Proportional Counter (LAXPC) spectrum. We explore various mechanisms that could possibly explain the presence of QPOs exclusively during the spin-down state.

1 INTRODUCTION

X-ray pulsars are highly magnetized neutron stars in binary systems, accreting matter from a companion star. Most X-ray pulsars are found in high-mass X-ray binaries (HMXBs) and possess strong magnetic fields of $\gtrsim {10}^{11}$ G. As the infalling material is channelled along the neutron star’s magnetic field lines, it impacts the magnetic poles, producing periodic X-ray pulsations. A small subset of X-ray pulsars resides in low-mass X-ray binaries (LMXBs), where the companion has a mass of $\lesssim 2{\mathrm{M}}_{\odot }$⁠. Only five such systems are currently known: Her X−1, 4U 1626−67, GX 1+4, GRO J1744−28, and 4U 1822−37 (see Bildsten et al. 1997; Jonker & van der Klis 2001). In contrast, the majority of accreting neutron stars in LMXBs have much weaker magnetic fields (⁠$B\lesssim {10}^9$ G), insufficient to strongly influence accretion, and are either non-pulsating or observed as accreting millisecond X-ray pulsars (AMXPs). Unlike AMXPs, which exhibit spin periods of a few milliseconds, X-ray pulsars typically have much longer spin periods, ranging from fractions of a second to thousands of seconds.

4U 1626−67 is an X-ray pulsar in an ultra-compact LMXB system, which was discovered in 1972 (Giacconi et al. 1972), with a short binary orbit of $\sim$42 min and a spin period of $\sim$7.7 s (Rappaport et al. 1977; Middleditch et al. 1981; Chakrabarty 1998; Jain et al. 2007). An ultra-compact environment in the presence of a strong magnetic field makes 4U 1626−67 an interesting source for studying accretion dynamics. The most prominent characteristic of 4U 1626−67 is its steady spin-up and spin-down episodes on time-scales up to several years (Chakrabarty et al. 1997; Jain & Paul 2009; Camero-Arranz et al. 2010; Jain, Paul & Dutta 2010; Sharma, Jain & Paul 2023b). During 1972–1990, 4U 1626−67 was steadily spinning up (labelled as 1SU, following the nomenclature of Sharma et al. 2023b). A torque reversal to spin-down (labelled as 1SD, hereafter) occurred around 1990 (Chakrabarty et al. 1997) and lasted around 18 yr. In 2008, 4U 1626−67 again showed a torque reversal to spin-up (labelled as 2SU, hereafter; Jain & Paul 2009; Camero-Arranz et al. 2010; Jain et al. 2010). In 2023, the source showed another torque reversal to spin-down (labelled as 2SD; Jenke, Wilson-Hodge & Malacaria 2023; Sharma et al. 2023b).

Since its discovery, coherent X-ray pulsations have always been detected in 4U 1626−67. However, all three torque reversals in this source have been accompanied by significant changes in the timing and spectral signatures at X-ray energies. The pulse profile has been strongly dependent on the torque state and energy (Jain et al. 2010; Beri et al. 2014; Sharma et al. 2023b). This is indicative of changes in the anisotropic radiative transfer within the accretion column in the presence of a strong magnetic field of $\sim {10}^{12}-{10}^{13}$ G (Kii et al. 1986). Over the years, the power density spectrum (PDS) has also displayed torque-state dependent occurrence of mHz quasi-periodic oscillations (QPOs, Kaur et al. 2008; Jain et al. 2010; Beri et al. 2014; Sharma et al. 2023b). Using Ginga data during 1SU, Shinoda et al. (1990) reported weak 0.04 Hz QPOs. During 1SD, strong $\sim$48 mHz QPOs were reported from observations made with instruments on-board Beppo-SAX (Owens, Oosterbroek & Parmar 1997), RXTE (Chakrabarty 1998; Kommers, Chakrabarty & Lewin 1998) and XMM–Newton (Krauss et al. 2007). During the spin-up phase, the QPO centroid frequency evolved with time from about 36 mHz in 1983 to about 49 mHz in 1993. But during the subsequent spin-down phase, it decreased at a rate of about 0.2 mHz yr${}^{-1}$ (Kaur et al. 2008). QPOs were not detected in the PDS during 2SU (Jain et al. 2010), however, XMM–Newton observation revealed both a 3 and a 48 mHz QPO, though the latter had a low rms amplitude (Beri, Paul & Dewangan 2018). In the current spin-down episode, $\sim$46 mHz QPOs re-appeared as reported from works based on the Nuclear Spectroscopic Telescope ARray (NuSTAR) observations (Sharma et al. 2023b; Tobrej et al.2024).

The torque reversals in 4U 1626−67 have also been associated with changes in the X-ray intensity. During 1SU, the source luminosity was about 10${}^{37}$ erg s${}^{-1}$ (White, Swank & Holt 1983). The X-ray flux decreased during the transition to 1SD (Chakrabarty et al. 1997), and the source luminosity increased by a factor of about three during the transition to 2SU (Jain et al. 2010). During the current phase (2SD), source luminosity dropped by a factor of $\sim$3 to the luminosities of $\sim 4\times {10}^{36}$ erg s${}^{-1}$ (Sharma et al. 2023b; Tobrej et al. 2024).

The X-ray continuum of 4U 1626−67 also correlates with the torque state. The spectrum during 1SU was described with a model comprising of $\sim$0.6 keV blackbody component and power-law component having photon index of $\sim$1 (Pravdo et al. 1979; Kii et al. 1986; Angelini et al. 1995). The source exhibited a similar spectral shape during the second spin-up phase (Jain et al. 2010; Camero-Arranz et al. 2012; Koliopanos & Gilfanov 2016; Sharma et al. 2023b). The energy spectrum was relatively hard in between these states (1SD, Owens et al. 1997; Vaughan & Kitamoto 1997; Yi & Vishniac 1999). The spectrum also shows the presence of strong emission lines of highly ionized Ne and O. The iron line is observed to be both torque-state and pulse-phase dependent (Beri, Paul & Dewangan 2015; Koliopanos & Gilfanov 2016; Beri et al. 2018). During the current 2SD, the broad-band spectrum is described by the empirical negative and positive power law with exponential cut-off (NPEX) model along with a soft blackbody temperature of $\sim$0.25 keV and the presence of cyclotron resonance scattering feature (CRSF) at $\sim$36 keV (Sharma et al. 2023b; Tobrej et al. 2024).

In the current work, we present the results from the timing and spectral analysis of an AstroSat target-of-opportunity (ToO) observation of 4U 1626−67 conducted during the 2SD torque state. The results are compared with the previous AstroSat observation of the source made during the 2SU (see Table 1). Additionally, we incorporate NuSTAR observations for timing analysis. The paper is organized in the following way. Section 2 provides a description of the observation and the data reduction procedure. In Section 3, the results from the timing and spectral analysis of 4U 1626−67 are presented. Our findings are discussed in Section 4.

2 OBSERVATION AND DATA REDUCTION

AstroSat is India’s first multiwavelength astronomical mission, which was launched in September 2015 by the Indian Space Research Organization (Agrawal 2006). It comprises of five scientific instruments, four of which can observe a source simultaneously from optical to hard X-ray energy bands (Agrawal 2006; Singh et al. 2014).

The timing and spectral analysis in this work has been performed using data from two observations of 4U 1626−67 taken with the Large Area X-ray Proportional Counter (LAXPC) onboard AstroSat (see Table 1). The first observation (obsID: 9000002100) was conducted on 2018 May 15 (MJD 58253), during the source’s spin-up phase (2SU), with a total exposure of 81 ks. The second observation (obsID: 9000005642) was performed on 2023 May 18 (MJD 60082), after the torque reversal to the spin-down phase (2SD), with an exposure of 45 ks. LAXPC has a timing resolution of 10 $\mu$s, and it consists of three co-aligned proportional counters (LAXPC10, LAXPC20, and LAXPC30) covering the energy range of 3–80 keV and has a total effective area of 6000 cm${}^2$ at 15 keV (Yadav et al. 2016; Agrawal et al. 2017). In the current analysis, we have used event analysis (EA) mode data from all layers of LAXPC20 only, as LAXPC30 was switched-off and LAXPC10 had variable gain due to gas leakage. The LAXPC software¹LaxpcSoft: version 3.4.4 was used to process the level 1 data. The source and background light curves and spectra were extracted by using the tool laxpcl1. The backshiftv3 tool was used to correct the background products for the gain shift. The 2023 data required an additional background down-scaling of 2.4 per cent to match the source and background data above 50 keV. The light curves were extracted with a binsize of 0.1 s. Solar system barycentre correction was performed by using as1bary  ² tool by taking JPL DE405 ephemeris and source position RA (J2000) = 16${}^{\mathrm{h}}$ 32${}^{\mathrm{m}}$ 16${}^{\mathrm{s}}$⁠.79 and Dec. (J2000) = $-{67}^{\mathrm{o}}$ 27${}^{\prime }$ 39${}^{\prime \prime }$⁠.3 (Lin, Webb & Barret 2012).

4U 1626−67 was also observed with the Soft X-ray Telescope (SXT) during both the AstroSat observations in the photon counting mode. But, owing to a low duty cycle (⁠$\sim$20 per cent) and a time resolution of $\sim$2.3 s, we have not used SXT data in this work.

3 RESULTS

3.1 Timing analysis

We extracted 0.1 s light curves in the 3–30 keV energy range, optimized for maximum signal-to-noise ratio and searched for X-ray pulsations using the epoch folding method (Leahy et al. 1983). The spin period determined from both the AstroSat observations is mentioned in Table 1. The bootstrap method described in Boldin, Tsygankov & Lutovinov (2013) was used to estimate the error in the spin period by simulating 1000 light curves (Sharma et al. 2023a). These results are consistent with Fermi/Gamma-ray Burst pulsar Monitor³ (Malacaria et al. 2020).

The 3–30 keV light curves from both AstroSat observations were folded with their corresponding pulse period to obtain pulse profiles (Fig. 1). The energy-resolved pulse profiles were extracted in the 3–5, 5–7, 7–10, 10–15, 15–20, 20–30, and 30–45 keV energy ranges, aligned by their minima and overlaid in Fig. 1 for comparison. At low energies (up to 20 keV), the pulse morphology during 2023 is quite different from that during 2018. At energies up to 5 keV, the 2018 profile shows two peaks separated by a broad, shallow dip covering about 0.5 pulse phase. Then, up to about 10 keV, two shallow dips appear at the pulse phase of 0.4 and 0.7. Thereafter, the profile tends to become sinusoidal. This is consistent with the profiles of the spin-up phase reported by Beri et al. (2014) using RXTE data and by Iwakiri et al. (2019) from Suzaku and NuSTAR data. The 2023 pulse profile shows several substructures up to about 7 keV. Beyond this, the shape is largely sinusoidal, and the width of the main dip increases with energy. This is somewhat similar to profiles seen by Krauss et al. (2007) during the first spin-down and by Sharma et al. (2023b) during the current phase. We have found that the 3–30 keV pulsed fraction for the 2023 profile is higher (⁠$57.3\pm 1.3$ per cent) than that of 2018 profiles (⁠$39.4\pm 0.4$ per cent).

Fig. 2 shows the 3–30 keV PDS of 4U 1626−67 generated from 2018 and 2023 data. The light curves were divided into 819.2 s segments to calculate the Fourier transform. The resulting PDSs were averaged and rebinned geometrically by a factor of 1.02. The PDS were rms normalized with white noise level subtracted to obtain the rms fractional variability of the time series. Both PDS exhibited sharp peaks at the spin frequency of the pulsar (⁠$\sim$130 mHz) with multiple harmonics and red noise. In addition to coherent pulsations, narrow excess in power at $\sim$46 and $\sim$83 mHz were observed in the 2023 PDS, indicating the presence of QPO. The 83 mHz QPO is not a harmonic of the 46 mHz and can be attributed to the beat between the 130 mHz coherent pulsation and the 46 mHz QPO (⁠${\nu }_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{t}}={\nu }_{\mathrm{N}\mathrm{S}}-{\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}$⁠). This is the first-ever detection of 83 mHz beat QPO during the current spin-down phase (2SD) of 4U 1626−67.

The PDS can be described with a sum of Lorentzian functions (e.g. Belloni, Psaltis & van der Klis 2002; Reig 2008; Sharma et al. 2024), each defined as

where ${\nu }_0$ is the centroid frequency, $\mathrm{\Delta }$ is the full-width at half-maximum (FWHM), and r is the integrated fractional rms The quality factor of Lorentzian $Q={\nu }_0/\mathrm{\Delta }$ is used to validate the presence of a QPO above the noise level. The Lorentzian components with $Q>2$ are generally considered as QPOs, which are otherwise labelled as band-limited noise (e.g. Belloni et al. 2002). The PDS of 4U 1626−67 from both observations in the 0.001–1 Hz range can be well described with two Lorentzian components representing the red noise continuum. For the 2023 observation, two additional Lorentzians were required to model QPOs at $\sim$46 and $\sim$83 mHz. Modelling of these features yielded a quality factor of 5.4 and 8.8, and rms values of $\sim$20 and $\sim$8 per cent, respectively. During the fitting, data points around the fundamental spin frequency and its harmonic were excluded. The 83 mHz QPO was detected with a 3.6 $\sigma$ significance, calculated by dividing the normalization of the Lorentzian component by its 1$\sigma$ lower bound. We also confirmed a false alarm probability of $<{10}^{-4}$ using simftest, corresponding to $>$3$\sigma$ significance. Table 2 summarizes the best-fitting parameters of PDS for both observations. Notably, the red noise component showed significant differences between the two epochs: the 2018 observation exhibited higher power at low frequencies and a broader pulse frequency peak compared to 2023. The low-frequency noise component had an rms of $\sim$26 per cent in 2018 and $\sim$15 per cent in 2023, while the high-frequency noise component showed an inverse trend, with $\sim$15 per cent rms in 2018 increasing to $\sim$38 per cent in 2023. This suggests a notable evolution in the source’s timing properties between the two epochs.

Fig. 3 shows the variation of the rms amplitude of 46 mHz QPO and 83 mHz beat QPO with energy for the 2023 observation. The horizontal error bar corresponds to the energy ranges used to generate the PDS, while the vertical error bars show the 1$\sigma$ error in the rms amplitude. The rms amplitude of fundamental 46 mHz QPO shows a clear positive correlation with energy. In contrast, the 83 mHz QPO shows no significant energy dependence.

Fig. 4 shows the residuals after modelling the red noise and 46 mHz QPO in the PDS for the 3–30 keV (top) and 15–20 keV (bottom) bands. The fundamental spin frequency and its harmonics are also indicated with their orders (green dashed lines) and the expected sidebands (red arrows) as if the QPO signal modulates the amplitude of the coherent pulsations. In addition to the 83 mHz feature, an upper sideband at 177$\pm$2 mHz (⁠${\nu }_{\mathrm{N}\mathrm{S}}+{\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}$⁠) was also apparent in both energy ranges. Energy-resolved PDS analysis revealed that the 177 mHz feature was significant only in the 15–20 keV range, with a detection significance of more than $3\sigma$ and rms of 10$\pm$3 per cent. The 15–20 keV PDS also showed additional features at other sideband frequencies, though at lower significance (⁠$<3\sigma$⁠). For instance, a feature at 615(12) mHz (lower sideband of n= 5) was detected with $2.5\sigma$ significance and rms of 12$\pm$4 per cent. Beyond 1 Hz (⁠$n>7$⁠), most sidebands appeared at a low confidence level in both energy bands, with no statistically robust detections, except for a prominent feature at 1.504(4) Hz. This feature was observed only in the 3–30 keV PDS, with a $4\sigma$ significance and rms of 9$\pm$2 per cent, and was not apparent in any of the energy-resolved PDS. The right panel of Fig. 4 shows the residuals near this frequency. A similar QPO at 1.3 Hz was previously observed in HMXB pulsar XTE J0111.2−7317 (Kaur et al. 2007), although the underlying origins are likely different. The 1.5 Hz feature in 4U 1626−67 could result from an overlap between n=11 upper sideband and n=12 lower sideband or may arise from instabilities in the accretion flow (Li et al. 2024).

3.1.1. NuSTAR PDS

We also conducted a search for sidebands in the NuSTAR (Harrison et al. 2013) light curves to investigate their origin. The analysis focused on the five NuSTAR observations carried out during the 2SD phase (Table A1). Details of these observations and the corresponding data reduction processes are outlined in Appendix  A. Initial results from the spectral and timing analyses of the 2SD phase were reported by Sharma et al. (2023b). Subsequently, Tobrej et al. (2024) presented spectral and timing results from four of the NuSTAR observations taken during 2SD. However, these observations have not been previously utilized to identify sidebands in the PDS.

For comparison with AstroSat, we extracted PDS in the 3–30 keV (broad band), 3–15 (low energy) and 15–30 keV (high energy) ranges from NuSTAR observations. Fig. 5 presents the residuals after modelling the NuSTAR PDS by the red noise and 46 mHz QPO. In obs 1, a feature at 177 mHz was clearly detected with $3\sigma$ significance and rms of 4.6$\pm$1.3 per cent in the 3–30 keV energy range. This feature was also present in the 3–15 and 15–30 keV range with lower significance, exhibiting rms of ${4.2}_{-1.5}^{+1.0}$ per cent and ${8.6}_{-3.5}^{+2.4}$ per cent, respectively. Additionally, a lower sideband peaked around 77 mHz and was dominant in 3–15 keV with rms of 7$\pm$2 per cent. In obs 2, which was taken just one day before the AstroSat observation, the 83$\pm$2 mHz sideband was detected with rms of 7$\pm$3 per cent, consistent with AstroSat but at a lower significance of $2.5\sigma$⁠, and only in the 3–15 keV band. No other sidebands could be detected significantly.

In obs 3, no sidebands were significantly detected in the 3–30 and 3–15 keV ranges. However, a feature at 191 mHz was present, albeit with low significance. In the 15–30 keV PDS, the 177(3) mHz feature was detected with $>3\sigma$ significance and rms of 12$\pm$3 per cent, consistent with AstroSat PDS in the 15–20 keV range. Additionally, an excess noise component was observed at 19$\pm$2 mHz with rms of $9_{-3}^{+4}$ per cent, which does not correspond to a sideband of any harmonics. In obs 4 and 5, no sidebands were detected. However, a feature at 92 mHz was identified with rms of ${6.1}_{-1.6}^{+1.3}$ per cent and a significance of $3\sigma$ in obs 4. This frequency is slightly offset from the expected 83 mHz sideband, suggesting the presence of additional oscillations in the source intensity.

3.2 Spectral analysis

To quantify the spectral variation, we performed a spectral analysis using LAXPC data from the 2018 and 2023 observations. A systematic uncertainty of 1.5 per cent was applied during spectral fitting. For the 2023 observation, we used data up to 45 keV, as the spectrum above this energy was dominated by background noise, while the 2018 spectrum extended up to 80 keV. Both spectra were rebinned using grppha to have a minimum count of 20 counts per bin and fitted both spectra simultaneously with an absorbed NPEX model (Mihara, Makishima & Nagase 1995; Makishima et al. 1999), following the methodology of Sharma et al. (2023b). The line-of-sight absorption, modelled with tbabs, was fixed at $N_H=9.6\times {10}^{20}$ cm${}^{-2}$ (HI4PI Collaboration 2016).

The absorbed NPEX model failed to provide a satisfactory fit, showing clear residuals around 37 keV, due to CRSF (Orlandini et al. 1998; Coburn et al. 2002; Iwakiri et al. 2019; Sharma et al. 2023b). To account for this, we included the cyclabs model, which features a pseudo-Lorentzian optical-depth profile to describe the CRSF (Makishima et al. 1990; Mihara et al. 1990). For the 2023 spectra, the CRSF line energy and width cannot be reliably constrained and, hence, were fixed at the values obtained from the 2018 observation. Additionally, an Fe K$\alpha$ emission line at 6.75 keV with a fixed width of 0.15 keV was included, following Sharma et al. (2023b), as allowing the line energy to vary resulted in the parameter pegging at 7 keV. An extra edge component at 10 keV (Manikantan, Paul & Rana 2023) was also required to fully describe the continuum, yielding a satisfactory fit.

The best-fitting model with their respective residuals is shown in Fig. 6, and the best-fitting spectral parameters are presented in Table 3. No significant spectral variability was detected between the two spectra from the different torque states. Although the average unabsorbed flux in the 3–30 keV range dropped from $7.7\times {10}^{-10}$ erg cm${}^{-2}$ s${}^{-1}$ in 2018 to $2.4\times {10}^{-10}$ erg cm${}^{-2}$ s${}^{-1}$ in 2023. This factor of $\sim$3 reduction in flux is consistent with previous findings by Sharma et al. (2023b).

4 DISCUSSION AND CONCLUSIONS

In this work, we have presented the results from the AstroSat observation of the ultra-compact X-ray pulsar 4U 1626−67, made after the recent torque reversal phase (2SU to 2SD) along with a comparison to the previous AstroSat observation from the 2SU phase. We also report results from the timing analysis of NuSTAR observations from the 2SD phase. The pulse profile follows the previously known torque-state dependence patterns (Beri et al. 2014; Sharma et al. 2023b). The flux in the 2SD dropped by a factor of three compared to 2SU. The CRSF was detected at 37 keV, consistent with previous reports (Orlandini et al. 1998; Coburn et al. 2002; Iwakiri et al. 2019; Sharma et al. 2023b). 4U 1626−67 is one of the few sources that exhibit both a CRSF and QPOs (Raman et al. 2021). This detection adds to the growing list of CRSF observations made with LAXPC (e.g. Varun et al. 2019a,b; Raman et al. 2021; Devaraj et al. 2024).

The PDS also exhibited torque-state dependence, with notable differences in noise variability and the appearance of QPOs, which were exclusive to the spin-down state (see Jain et al. 2010; Sharma et al. 2023b). The base of the pulse peak in the PDS from 2SU (2018) shows a broadened shape compared to 2SD (2023), consistent with Jain et al. (2010). This behaviour is similar to X-ray pulsar 4U 1901+03, where the pulse peak narrows when the QPO feature begins to emerge (James et al. 2011). However, it contrasts with IGR 19294+1816, where the QPO feature appears alongside the broadened peak (Raman et al. 2021). This broadening of pulse peak arises from the coupling between the coherent periodic pulse variability and the aperiodic red noise component (Lazzati & Stella 1997). The increased power at low frequencies in the PDS during the spin-up state can be coupled with the pulse variability, leading to the observed broadening. Similar behaviour is observed in 4U 1901+03, where increased low-frequency variability coincides with pulse peak broadening (James et al. 2011). When the QPO appears, low-frequency noise is suppressed, resulting in less broadening of the pulse peak.

A pronounced $\sim$46.5 mHz QPO feature with $\sim$20 per cent rms and quality factor of $\sim$5 was detected during the 2023 observation (2SD). Additionally, a secondary QPO was observed at $\sim$83.6 mHz, having a quality factor of 8.8 and rms of 7.8 per cent. The appearance of QPOs in X-ray pulsars is often explained by two primary models: the Keplerian Frequency Model (KFM, van der Klis et al. 1987) and the Beat Frequency Model (BFM, Alpar & Shaham 1985; Lamb et al. 1985). As per the KFM, inhomogeneous structures in the Keplerian disc regularly attenuate the pulsar beam. In the case of 4U 1626−67, the option of KFM is ruled out because that requires the neutron star to spin faster than the inner accretion disc (⁠${\nu }_{\mathrm{N}\mathrm{S}}>{\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}$⁠), which will result in centrifugal inhibition of mass accretion (Stella, White & Rosner 1986). In the BFM, the QPO is attributed to the modulation in mass accretion rate onto the neutron star’s poles at the beat frequency between the spin frequency of the neutron star (⁠${\nu }_{\mathrm{N}\mathrm{S}}$⁠) and the Keplerian orbital frequency (⁠${\nu }_{\mathrm{k}}$⁠) of the inner accretion disc. According to BFM, ${\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}={\nu }_{\mathrm{k}}-{\nu }_{\mathrm{N}\mathrm{S}}$⁠. The 46 mHz QPO in 4U 1626−67 exhibits a torque-state dependence, apparent during the spin-down phase but absent during the spin-up phase. A weak 0.04 Hz QPO was reported during 1SU (Shinoda et al. 1990), but it was never detected during 2SU. Therefore, its reappearance during the recent torque reversal to the spin-down state clearly establishes its linkage to the torque state. Since ${\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}<{\nu }_{\mathrm{N}\mathrm{S}}$⁠, the BFM can provide a plausible explanation for the 46 mHz QPO, corresponding to Keplerian frequency of $\sim$177 mHz and inner disc radius of $5.3\times {10}^8$ cm (Manikantan et al. 2024).

However, the detection of an additional QPO at 83 mHz during the current 2SD phase with AstroSat adds a layer of complexity to the BFM. This frequency is not a harmonic of the primary 46 mHz QPO. Interestingly, the 83 mHz QPO follows the relationship ${\nu }_{83\mathrm{m}\mathrm{H}\mathrm{z}}+{\nu }_{46\mathrm{m}\mathrm{H}\mathrm{z}}={\nu }_{\mathrm{N}\mathrm{S}}$⁠, indicative of a beat between the (46 or 83 mHz) QPO and spin of neutron star, or the lower sideband (⁠${\nu }_{\mathrm{l}\mathrm{s}}={\nu }_{\mathrm{N}\mathrm{S}}-{\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}$⁠) of the spin frequency. Similar sidebands were observed in the 1SD phase (Kommers et al. 1998), caused by the modulation of coherent pulsation by the amplitude of the QPOs. Kommers et al. (1998) suggested that these QPOs are produced by a structure orbiting the neutron star at the QPO frequency rather than originating at the inner edge of the accretion disc. However, it is not clear how such a structure could persist in the accretion disc against the differential rotation for more than a decade (during 1SD), disappear during the 2SU phase, and then reappear in the current 2SD phase. If we consider the QPO frequency to represent the inner disc Keplerian frequency ${\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}={\nu }_{\mathrm{k}}$⁠, it implies that the accretion disc lies outside the corotation radius (⁠$R_{\mathrm{c}\mathrm{o}}$⁠) of the neutron star. In this scenario, the QPOs could be generated by the interaction between the magnetic field of the neutron star and the inner edge of the accretion disc, where a portion of the disc infiltrates the corotation radius due to a misalignment of the magnetic axis. Notably, the X-ray flux of 4U 1626−67 has been observed to decrease by a factor of $\sim$3 during the current spin-down phase (see Sharma et al. 2023b), further strengthening the hypothesis of the disc’s location outside or close to the corotation radius.

Other models, like thermal disc instabilities and magnetic disc precession models, are unlikely to account for the observed frequencies, as they typically predict variability of the order of $\sim$1 mHz (Shirakawa & Lai 2002; Roy et al. 2019; Liu et al. 2022). Another promising explanation is the model proposed by D’Angelo & Spruit (2010), which has been used to describe the 40 mHz QPO in Cen X-3 (Liu et al. 2022). In this scenario, the disc is truncated just outside $R_{\mathrm{c}\mathrm{o}}$⁠, and centrifugal forces inhibit accretion by transferring angular momentum from the neutron star to the disc, causing the neutron star to spin-down. Instead of the mass being ejected, the accreting mass can stay piled up at high surface density in the inner disc, just outside corotation. Over time, mass accumulates in the disc, increasing gas pressure, which eventually pushes the inner edge of the disc inside $R_{\mathrm{c}\mathrm{o}}$⁠, allowing accretion to resume. Once accretion occurs, the gas pressure decreases, and the disc moves back outside $R_{\mathrm{c}\mathrm{o}}$⁠, restarting the cycle. The period of this cycle can range from 0.02 to 20 viscous time-scales (⁠$t_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}$⁠) of a few hundred seconds (D’Angelo & Spruit 2010). This mechanism could provide a reasonable explanation for the observed mHz QPOs in 4U 1626−67 exclusively to the spin-down torque state.

The 46 mHz QPO exhibited strong energy dependence, with its rms increasing with energy (see Manikantan et al. 2024), while the 83 mHz QPO did not show any significant energy dependence (Fig. 3). However, the upper sideband at 177 mHz displayed clear energy dependence, being detected exclusively in the 15–20 keV energy range, whereas the lower sideband at 83 mHz was observed across all selected energy bands with AstroSat. The non-detection of the 177 mHz feature above 20 keV is likely attributed to the increased background noise in LAXPC at higher energies. The NuSTAR observations further confirmed the energy-dependent nature of these sidebands, with the rms of the 177 mHz sideband increasing from $\sim$4  per cent in the 3–15 keV range to $\sim$8.6  per cent in the 15–30 keV range during obs 1. In obs 3, the 177 mHz sideband was detected only in the 15–30 keV range with an rms of $\sim$12  per cent. These sidebands also showed time dependence, which could be due to the statistical limitations of short observation span. During the 1SD phase, Kommers et al. (1998) reported both lower and upper sidebands. These sidebands were symmetric in frequency, but their power amplitude was energy-dependent. These sidebands arise due to the Fourier frequency-shifting theorem, where the QPO signal modulates the amplitude of the coherent pulsations, producing symmetric sidebands in the shape of the QPO peak around the harmonic frequencies $n{\nu }_{\mathrm{N}\mathrm{S}}\pm {\nu }_{\mathrm{Q}\mathrm{P}\mathrm{O}}$⁠, where n=1, 2, … is an integer (Kommers et al. 1998). If the sidebands in 4U 1626−67 are due to amplitude modulation of the 46 mHz QPO, there must be an additional mechanism that can explain the enhancement or reduction of the power amplitude of these sidebands in different energy ranges. Interestingly, the lower sideband is almost absent in weakly magnetized neutron star binaries 4U 1608–52 and 4U 1728–34 (Jonker, Méndez & van der Klis 2000). In these cases, the magnetic field influence on the accretion process is weaker, suggesting that the interaction between the neutron star’s magnetic field and the accretion disc in 4U 1626−67 plays a crucial role in generating the sidebands and preferentially amplifying one sideband.

ACKNOWLEDGEMENTS

This work has made use of data from the AstroSat mission of the Indian Space Research Organisation (ISRO), archived at the Indian Space Science Data Centre (ISSDC). We thank the LAXPC Payload Operation Center (POC) for verifying and releasing the data via the ISSDC data archive and providing the necessary software tools. We have also made use of the software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This research has also made use of NuSTAR data obtained from the HEASARC and the NuSTAR Data Analysis Software (NUSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). CJ acknowledges the financial assistance received from the ANRF (erstwhile SERB) – DST grant (CRG/2023/000043). AB acknowledges SERB (SB/SRS/2022-23/124/PS) for financial support and is also grateful to the Royal Society, UK. We also thank the anonymous referee for insightful comments and suggestions.

DATA AVAILABILITY

Data used in this work can be accessed through the Indian Space Science Data Center (ISSDC) at https://astrobrowse.issdc.gov.in/astro_archive/archive/Home.jsp and the HEASARC archive at https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl.

Footnotes

REFERENCES

APPENDIX A: NuSTAR OBSERVATIONS AND DATA REDUCTION

4U 1626−67 was observed five times by NuSTAR (Harrison et al. 2013) following the torque reversal in 2023 (Table A1). The data were processed following the procedures outlined by Sharma et al. (2023b), utilizing heasoft version 6.34 along with the latest calibration files (version 20250122). The calibrated and screened event files were generated using the task nupipeline. A circular region of radius 80 arcsec centred at the source position was used to extract the source events. Background events were extracted from a circular region of the same size away from the source. The task nuproduct was used to generate the science products such as light curves. These light curves were corrected to the Solar system barycentre using barycorr. Finally, the background-corrected light curves from the focal plane modules FPMA and FPMB were combined using lcmath to obtain the summed light curves for subsequent timing analysis.