# Recovery of High-energy Low-frequency Quasi-periodic Oscillations from Black
Hole X-ray Binary MAXI J1535–571 with a Hilbert-Huang Transform Method

Qing C. Shui Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, 100049, Beijing, China  University of
Chinese Academy of Sciences, Chinese Academy of Sciences, 100049, Beijing,
China  S. Zhang  Key Laboratory of Particle Astrophysics, Institute of High
Energy Physics, Chinese Academy of Sciences, 100049, Beijing, China  Shuang N.
Zhang Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, 100049, Beijing, China  University of
Chinese Academy of Sciences, Chinese Academy of Sciences, 100049, Beijing,
China  Yu P. Chen  Key Laboratory of Particle Astrophysics, Institute of High
Energy Physics, Chinese Academy of Sciences, 100049, Beijing, China  Ling D.
Kong Institut für Astronomie und Astrophysik, Kepler Center for Astro and
Particle Physics, Eberhard Karls, Universität, Sand 1, D-72076 Tübingen,
Germany  Jing Q. Peng Key Laboratory of Particle Astrophysics, Institute of
High Energy Physics, Chinese Academy of Sciences, 100049, Beijing, China
University of Chinese Academy of Sciences, Chinese Academy of Sciences,
100049, Beijing, China  L. Ji  School of Physics and Astronomy, Sun Yat-Sen
University, Zhuhai, 519082, China  Peng J. Wang Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
100049, Beijing, China  University of Chinese Academy of Sciences, Chinese
Academy of Sciences, 100049, Beijing, China  Z. Chang Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 100049, Beijing, China  Zhuo L. Yu  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
100049, Beijing, China  Hong X. Yin Shandong Key Laboratory of Optical
Astronomy and Solar-Terrestrial Environment, School of Space Science and
Physics, Institute of Space Sciences, Shandong University, Weihai, Shandong
264209, China  Jin L. Qu Key Laboratory of Particle Astrophysics, Institute of
High Energy Physics, Chinese Academy of Sciences, 100049, Beijing, China  L.
Tao Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, 100049, Beijing, China  Ming Y. Ge Key Laboratory
of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 100049, Beijing, China  X. Ma Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
100049, Beijing, China  L. Zhang  Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 100049,
Beijing, China  W. Yu  Key Laboratory of Particle Astrophysics, Institute of
High Energy Physics, Chinese Academy of Sciences, 100049, Beijing, China  J.
Li  CAS Key Laboratory for Research in Galaxies and Cosmology, Department of
Astronomy, University of Science and Technology of China, Hefei 230026, China

###### Abstract

We propose a method based on the Hilbert-Huang transform (HHT) to recover the
high-energy waveform of low-frequency quasi-periodic oscillations (LFQPOs).
Based on the method, we successfully obtain the modulation of the phase-folded
light curve above 170 keV using the QPO phase reconstructed at lower energies
in MAXI J1535–571 with Insight-HXMT observations. A comprehensive simulation
study is conducted to demonstrate that such modulation indeed originates from
the QPO. Thus the highest energies turn out to significantly exceed the upper
limit of $\sim$100 keV for QPOs reported previously using the Fourier method,
marking the first opportunity to study QPO properties above 100 keV in this
source. Detailed analyses of these high-energy QPO profiles reveal different
QPO properties between the 30–100 keV and 100–200 keV energy ranges: the phase
lag remains relatively stable, and the amplitude slightly increases below
$\sim 100$ keV, whereas above this threshold, soft phase lags and a decrease
in amplitude are observed. Given the reports of a hard tail detection in broad
spectroscopy, we propose that the newly discovered QPO properties above 100
keV are dominated by the hard tail component, possibly stemming from a
relativistic jet. Our findings also indicate a strong correlation between the
QPOs originating from the jet and corona, supporting the scenario of jet-
corona coupling precssion. We emphasize that our proposed HHT-based method can
serve as an efficient manner in expanding the high energy band for studying
QPOs, thereby enhancing our understanding of their origin.

Accretion (14) — Black hole physics (1736) — X-ray binary stars (1811) —
Stellar mass black holes (1611)

##  1 Introduction

Following a prolonged period of quiescence at low-flux levels, black hole
X-ray binaries (BHXRBs) can experience outbursts, with substantial changes in
their X-ray emission properties (Remillard & McClintock, 2006; Done et al.,
2007). During a complete outburst, the source typically exhibits four distinct
canonical states: the low/hard state (LHS), hard intermediate state (HIMS),
soft intermediate state (SIMS) and high soft state (HSS). A source exhibits
distinct X-ray spectral and variability properties in different states
(Belloni et al., 2005; Homan & Belloni, 2005). Low frequency quasi-periodic
oscillations (LFQPOs) are observed in most transient BHXRBs, typically ranging
from 0.1 to 30 Hz (van der Klis, 1989). These oscillations are characterized
by a narrow peak with finite width in the power density spectra (PDS) and are
typically classified into three types (A, B, and C) based on their centroid
frequency, $Q$ factor (ratio of frequency to the full width at half maximum),
and root-mean-square (rms) amplitude (Wijnands et al., 1999; Casella et al.,
2005; Remillard & McClintock, 2006).

LFQPOs are thought to originate from the innermost part of the accretion flow,
where the most X-ray emissions and fast variability are produced. Fast X-ray
variability is observed across a wide range of timescales (Belloni & Stella,
2014). In addition to the QPO component, broadband continuum components also
contribute to the PDS, which represent fast aperiodic variability. Several
models have been suggested to account for this noise component, such as the
shot noise model, coronal flare model, and fluctuation propagation model
(Terrell, 1972; Nolan et al., 1981; Scargle et al., 1993; Mineshige et al.,
1994; Ingram & Done, 2011). For LFQPOs, various models have been also
proposed, based on either the geometric effect or the intrinsic variability of
the accretion flow. Intrinsic models include the trapped corrugation modes
(Kato, 1990; Wagoner, 1999), the Accretion-ejection instability model (AEI,
Tagger & Pellat, 1999), and the Two-Component Advection Flow model (Molteni et
al., 1996), etc. Most geometric models, on the other hand, are associated with
the relativistic Lense-Thirring (L-T) precession effect, such as the
relativistic precession model (RPM; Stella et al., 1999), L-T precession of
the entire hot accretion flow (Ingram et al., 2009) and precession of the jet
(Stevens & Uttley, 2016; de Ruiter et al., 2019; Ma et al., 2021, 2023).

Componization is widely considered to be the radiative mechanism behind type-C
QPOs (e.g., Méndez et al., 2022; Bellavita et al., 2022), because
observational studies show that the factional rms of type-C QPOs generally
increases with photon energy (Zhang et al., 2017; Huang et al., 2018; Kong et
al., 2020; Zhang et al., 2020a) and no evident disc-like component is detected
in their rms spectra (Sobolewska & Życki, 2006; Axelsson et al., 2014;
Axelsson & Done, 2016; Shui et al., 2023a). Furthermore, Ma et al. (2021)
reported the discovery of tpye-C QPOs above 200 keV, further indicating a
strong association with Comptonized emission. Additionally, inclination-
dependent amplitudes and time lags of type-C QPOs (see Motta et al., 2015; van
den Eijnden et al., 2017), and QPO phase-dependent reflection component
revealed with phase-resolved spectroscopy (see Ingram & van der Klis, 2015;
Ingram et al., 2016, 2017; Nathan et al., 2022; Shui et al., 2023a), may
suggest type-C QPOs originate from some geometrical effects. For recent
reviews of observations and theories of LFQPOs, the reader is referred to
Ingram & Motta (2019).

The study of LFQPOs in high energy bands ($>30$ keV) is essential to our
understanding of their origin, since the radiative mechanism and geometric
properties of LFQPOs could be further constrained (Huang et al., 2018; Kong et
al., 2020; Ma et al., 2021; Liu et al., 2021). Thanks to the wide energy range
of Insigt-HXMT, we can study high-energy LFQPOs above 100 keV in some BHXRBs
in details (see e.g. Ma et al., 2021). Furthermore, $\gamma$-ray QPOs have
been detected by Fermi-LAT in numerous bright blazars (see e.g. Sandrinelli et
al., 2014; Ackermann et al., 2015; Banerjee et al., 2023). Consequently,
certain BHXRBs may have the potential to exhibit $\gamma$-ray QPOs, which
could involve different radiative mechanisms from those observed in the X-ray
bands, and warrant further detailed investigation in the future. In this
Letter, we propose a novel technique based on the Hilbert-Huang transform
(HHT) for recovering high-energy QPO signals in MAXI J1535–571, which is one
of the brightest BHXRBs observed by Insight-HXMT. The HHT is a powerful tool
for analyzing phenomena with non-stationary periodicity proposed by Huang et
al. (1998) and has been successfully applied in astronomical researches (Camp
et al., 2007; Hu et al., 2011, 2014; Su et al., 2015; Hsieh & Chou, 2020; Yu
et al., 2023a; Shui et al., 2023a).The HHT allows us to track changes in the
frequency of QPOs on short time-scales and conduct their phase-resolved
spectral analysis.

MAXI J1535–571 is an X-ray transient discovered by Gas Slit Camera of the
Monitor of All-sky X-ray Image (MAXI/GSC) and and the Burst Alert Telescope of
Swift (Swift/BAT) independently when it went into outburst on 2017 September 2
(Kennea et al., 2017; Negoro et al., 2017). Subsequently, this outburst was
extensively monitored at radio (Parikh et al., 2019; Russell et al., 2019),
sub-millimetre (Tetarenko et al., 2017), infrared and optical (Baglio et al.,
2018), and X-ray (Huang et al., 2018; Nakahira et al., 2018; Stevens et al.,
2018; Tao et al., 2018; Bhargava et al., 2019) wavelengths. From HI
absorption, the source distance was estimated to be $4.1^{+0.6}_{-0.5}$kpc
(Chauhan et al., 2019). LFQPOs (type-A, B and C) were detected by multiple
X-ray instruments in the spectral state of LHS, HIMS and SIMS (Huang et al.,
2018; Stiele & Kong, 2018; Mereminskiy et al., 2018; Stevens et al., 2018).

In this Letter, we propose a HHT-based technique to reconstruct the high-
energy QPO waveforms. In Section 2, we outline the Insight-HXMT observations
and data reductions. Then, in Section 3, we present a brief introduction to
the HHT method and its application in recovering high-energy QPO signals, and
then present the corresponding results derived from MAXI J1535–571. Finally,
we discuss and summarize these results in Section 4 and Section 5,
respectively.

##  2 Observations and Data Reduction

Launched on June 15, 2017, the Hard X-ray Modulation Telescope (Insight-HXMT)
is the first Chinese X-ray astronomy satellite (Zhang et al., 2014, 2020b).
The science payload of Insight-HXMT allows for observations across a broad
energy band (1–250 keV) using three telescopes: the High Energy X-ray
telescope (HE, composed of NaI/CsI, 20–250 keV), the Medium Energy X-ray
telescope (ME, with a Si pin detector, 5–30 keV), and the Low Energy X-ray
telescope (LE, using an SCD detector, 0.7–13 keV). For additional information
about the Insight-HXMT mission, please refer to Zhang et al. (2019); Liu et
al. (2020); Cao et al. (2020); Chen et al. (2020).

We focus on high-energy signals of the type-C QPO detected by Insight-HXMT, so
one of the bright systems MAXI J1535–571 with numerous type-C QPO observations
(see Huang et al., 2018; Kong et al., 2020; Zhang et al., 2022; Yu et al.,
2023b) is investigated in this study. We extract the event data using Insight-
HXMT Data Analysis Software v2.05111http://hxmtweb.ihep.ac.cn/SoftDoc.jhtml
and the current calibration model
v2.06222http://hxmtweb.ihep.ac.cn/software.jhtml, following the criteria
recommended by the Insight-HXMT team: (1) elevation angle (ELV) greater than
$10^{\circ}$; (2) geometric cutoff rigidity (COR) greater than 8 GeV; (3)
pointing position offset less than $0.04^{\circ}$; and (4) good time intervals
(GTIs) at least 300 s away from the South Atlantic Anomaly (SAA). Backgrounds
are derived from blind detectors using the LEBKGMAP, MEBKGMAP, and HEBKGMAP
tools, version 2.0.9, based on the standard Insight-HXMT background models
(Liao et al., 2020a, b; Guo et al., 2020). All photons are corrected to the
barycenter of the solar system for their arrival times using the HXMTDAS tool
hxbary. In Figure 1, we present the long-term net light curve for MAXI
J1535–571 obtained from HE in the energy range of 30–120 keV. To ensure robust
count statistics for HE, we select QPO observations with HE count rates in the
energy range of 30–120 keV that closely approximate the highest level across
the outburst (as indicated by the shaded area in Figure 1). We subsequently
generate net light curves with a time resolution of 0.05 s for these selected
observations, which are then subjected to the detailed timing analysis.
Detailed information regarding the aforementioned observations is available in
Table 1.

Figure 1: HE light curves of MAXI J1535–571. In this plot, red triangles denote observations with type-C QPO detection, while blue circles represent those without type-C QPO detection. Each point represents one observation, and the shaded area indicates the selected observations for detailed timing analyses. Table 1: Log of Insight-HXMT Observations Used in This Work # | ObsID | Start Time | QPO Frequency
---|---|---|---
|  | (MJD) | (Hz)
1 | P011453500144 | 58008.48 | $2.53\pm 0.02$
2 | P011453500145 | 58009.29 | $2.61\pm 0.01$
3 | P011453500201 | 58010.32 | $1.78\pm 0.01$
4 | P011453500301 | 58011.25 | $2.02\pm 0.01$
5 | P011453500401 | 58012.38 | $2.76\pm 0.01$
6 | P011453500501 | 58013.37 | $3.33\pm 0.01$
7 | P011453500601 | 58014.14 | $3.32\pm 0.03$

##  3 Analysis and Results

The Fourier transform is a general tool for searching QPO signals in X-ray
binaries. QPOs can be identified as a series of harmonic peaks with narrow
widths in the PDS. Due to limited count statistics in the high-energy bands
($>100$ keV), previous studies typically combined multiple observations to
achieve a higher signal-to-noise ratio (SNR) (see e.g., Ma et al., 2021; Yang
et al., 2022). Conversely, periodic folding is a common approach used in
studies of neutron star pulsations to search for high-energy pulse profiles
(see e.g., Hou et al., 2022). Since it has been demonstrated in BHXRBs that
there exists an average underlying waveform of the LFQPO (Ingram & van der
Klis, 2015; de Ruiter et al., 2019), similar folding methods could potentially
reconstruct the high-energy QPO waveform if the high-energy QPO signal indeed
exists. However, when dealing with QPO signals in black hole binaries, a
simple folding of the light curve on a period is not suitable due to the non-
deterministic evolution of their phases over time (see e.g., Ingram & van der
Klis, 2015). Specifically, Morgan et al. (1997) observed the QPO phase to
drift on a random walk away from that of a strictly periodic signal. There are
several well-known methods for studying the time variations of the QPO
frequency and amplitude, such as dynamic power spectra (DPS, see e.g. Homan et
al., 2020) and wavelet analysis (see e.g. Lachowicz & Done, 2010; Zhang et
al., 2023). In general, “instantaneous” frequencies obtained from the DPS or
wavelet analysis can be used for segmentally periodic folding the high-energy
light curve to construct high-energy QPO waveforms. However, the DPS requires
setting parameters for the width of the time window, and the wavelet analysis
needs parameters for the decay length of the wavelet function, both of which
significantly affect the performance of the methods. In addition to periodic
folding, phase-folding is another approach for constructing the waveform for a
non-stationary periodicity (see e.g. Su et al., 2015; Shui et al., 2023a).
Time-resolved complex power spectra can provide time variation of the phase
information for the QPO, but it also faces challenges related to the selection
of the width of the time window. Furthermore, it is still under debate how to
separate the intrinsic QPO phase from the complex Fourier transform, as both
QPO and noise components contribute to the Fourier transform at the QPO
frequency (Ma et al., 2021; Zhou et al., 2022). The HHT is a powerful timing
analysis tool that can provide an instantaneous phase function for the
intrinsic QPO component at each time bin (see e.g. Su et al., 2015; Shui et
al., 2023a), effectively addressing the challenges of the aforementioned
methods. In this section, we will first introduce a method for searching high-
energy QPOs based on the HHT method.

###  3.1 Hilbert-Huang Transform Analysis

The HHT, proposed by Huang et al. (1998) as an adaptive data analysis method,
offers a powerful tool for analyzing signals with non-stationary periodicity.
It consists of two main components: mode decomposition, which aims to
decomposing a signal into several intrinsic mode functions (IMFs), and Hilbert
spectral analysis (HSA), which allows for the extraction of frequency and
phase functions of the desired IMFs, such as the QPO component (e.g. Hu et
al., 2014; Hsieh & Chou, 2020). Following Shui et al. (2023a), we use the
variational mode decomposition (VMD, Dragomiretskiy & Zosso, 2014) method to
extract the IMF correspond to the QPO component. Compared to the traditional
method empirical mode decomposition (EMD, Huang et al., 1998; Huang & Wu,
2008) technique, the VMD method theoretically mitigates the influence of mode
mixing issue by decomposing the signal into a sum of IMFs with analytically
calculated limited center frequency and bandwidth. Mode mixing is one of the
major drawbacks of EMD, which is defined as a single IMF either consisting of
signals of widely disparate scales or a signal of a similar scale residing in
different IMF components (Huang & Wu, 2008). In the VMD algorithm, the set of
IMFs of an input signal, $f(t)$, is initialized as
$\\{u_{k}\\}=\\{u_{1},...,u_{K}\\}$ along with their corresponding center
frequencies $\\{\omega_{k}\\}=\\{\omega_{1},...,\omega_{K}\\}$. The objective
of the VMD is to minimize the sum of bandwidths of IMFs by solving the
following optimization problem333In Equation (1), $||h(x)||_{2}$ means the
$L^{2}$ (Euclidean) norm of the function $h(x)$, which is calculated as
$\left[\int{|h(x)|^{2}dx}\right]^{1/2}$. For a discrete vector $h_{k}(x_{k})$,
the $L^{2}$ norm of it is calculated as
$\left[\sum_{k=1}^{K}|h_{k}|^{2}\right]^{1/2}$. From another perspective, the
inner product can be expressed as $\langle
p(x),q(x)\rangle=\int{p^{*}(x)q(x)dt}$, which implies that
$||h(x)||_{2}=\left[\langle h(x),h(x)\rangle\right]^{1/2}$.:

$\min_{\\{u_{k}\\},\\{\omega_{k}\\}}\left\\{\alpha\sum_{k=1}^{K}{\rm
BW}_{k}^{2}+\left|\left|f(t)-\sum_{k=1}^{K}u_{k}(t)\right|\right|_{2}^{2}\right\\},$
(1)

where $K$ is the total number of modes, BWk is the bandwidth of the $k$th
mode. The first component of the objective function enhances mode compactness,
whereas the second component quantifies the reconstruction error, which ought
to be minimized. The parameter $\alpha$ acts as a balancing weight between
compactness and the reconstruction error. In Appendix A, we offer a concise
overview of the VMD algorithm. For the complete algorithmic description of the
VMD, we recommend referring to Dragomiretskiy & Zosso (2014). There are two
parameters should be initially set in the VMD algorithm, which is the total
number of modes $K$ and the weighting factor $\alpha$. In the VMD algorithm,
selecting the parameters $K$ and $\alpha$ is vital for its efficacy.
Generally, an increase in $K$ correlates with a rise in $\alpha$, suggesting
that the bandwidth narrows as the number of modes grows. Following the
procedure described in Shui et al. (2023a), we initially assign $\alpha$ a
typical value of 2000 and opt for 2 modes ($K=2$). We then generate the power
spectrum for each mode and compare it with the power spectrum of the original
light curve. Should no mode display a power spectrum mirroring the QPO
distinctive pattern, we increment $K$. This iterative strategy persist until
the fundamental and harmonic components of the QPO are discernibly isolated in
individual modes. Finally, we adjust $\alpha$ to align the IMF PDS profile
with the QPO peak in the power spectrum of the original light curve. The
selection of $K$ and $\alpha$ involves elements of subjectivity, which is one
of the major concerns of the VMD algorithm (Nazari & Sakhaei, 2018). This
issue requires future investigations, which are beyond the scope of this
study. Nevertheless, we find that making moderate adjustments to the values of
$K$ and $\alpha$ have a minimal impact on the subsequent Hilbert analysis when
ensuring that the PDS of the desired IMF aligns with the QPO peak in the
observed PDS (see Figure 2). In addition, the current VMD algorithm addresses
certain limitations of the original EMD algorithm, and show attractive
performance compared to the existing mode decomposition methods. The code we
used to perform VMD is from vmdpy v0.2444https://github.com/vrcarva/vmdpy
(Carvalho et al., 2020), which is an open-source Python package. It is worth
noting that an appropriate time resolution is required for the HHT analysis,
as too dense sampling results in poor SNR in each bin, while too sparse
sampling does not sufficiently decompose the QPO component. In our case, where
the QPO frequency is $\sim$2–3 Hz and the net count rate is $\sim 600$ cts
s-1, a time resolution of 0.05s is appropriate. This ensures that the time
internal between the samplings is much shorter than the variability timescale,
and each time bin has a signal-to-noise ratio (SNR) $\gtrsim 5\sigma$. Figure
2a shows an example of a 14-s-long light curve with a 0.05s time resolution
from HE observations of MAXI J1535–571 (obsID P011453500201) and the
corresponding IMFs. In this representative observation, the observed light
curve exhibits an LFQPO around 2 Hz. For the VMD algorithm, we selected
parameters with $K=5$ and $\alpha=300$. We successfully identify a $\sim$2 Hz
oscillation and its harmonic component as the second and third IMFs,
respectively. To compare results of the VMD algorithm with the Fourier
analysis, Figure 2b plots the power spectra for the original HE light curve
alongside that of the second and third IMFs. These power spectra are generated
in the frequency range of 0.04–10 Hz. A simulation also demonstrates the
robustness of the VMD algorithm that it shows good performance in the presence
of noise in the analyzed observations (see Appendix B for details).

Figure 2: Hilbert-Huang transform analysis of a representative example of a
14-s-long light curve with a 0.05s time resolution of HE in energy range of
30–120 keV. a, representative example of a 14-s-long light curve and
corresponding five IMFs obtained from VMD algorithm. b, power density spectra
of the original HE light curve, as well as the second and third IMFs. c, the
QPO (the second) IMF and corresponding phase function, instantaneous frequency
and amplitude. Figure 3: The Hilbert-Huang transform (HHT) phase-folded QPO
waveforms of MAXI J1535–571 from Insight-HXMT/HE light curves. To enhance the
SNR, waveforms from all selected observations have been collectively combined,
and the final combined waveform are normalized by the mean value of the count
rate. These QPO phases are obtained through HHT analysis of light curves in
the 30–120 keV energy band. Each panel showcases the significance level
determined via the cross-correlation technique, as detailed in Appendix C.

For each IMF obtained from the mode decomposition, the Hilbert transform can
calculates its physically meaningful instantaneous phase, amplitude and
frequency functions. For a time series $f(t)$, its Hilbert transform is
defined as

$\mathcal{H}\left[f(t)\right]=\frac{1}{\pi}{\rm
p.v.}\int{\frac{f(\tau)}{t-\tau}d\tau},$ (2)

where $\rm p.v.$ is the Cauchy principal value. Using the Hilbert transform,
one can obtain the corresponding analytic signal, which is defined as

$f_{\rm A}(t)=f(t)+j\mathcal{H}\left[f(t)\right]=A(t)e^{j\phi(t)},$ (3)

where $j$ is the imaginary unit. The time-dependent amplitude, $A(t)$, and
phase, $\phi(t)$, can be obtained from

$A(t)=\left\\{f^{2}(t)+\mathcal{H}^{2}\left[f(t)\right]\right\\}^{1/2}$ (4)

and

$\phi(t)=\arctan{\left\\{\frac{\mathcal{H}(\left[f(t)\right]}{f(t)}\right\\}},$
(5)

respectively. Therefore, the instantaneous frequency, $\nu(t)$, can be defined
as

$\nu(t)=\frac{1}{2\pi}\frac{d\phi(t)}{dt}.$ (6)

Considering the focus of the present study on investigating QPOs, we only
apply the Hilbert transform to the IMF corresponding to the QPO component
(e.g. the second IMF plotted in Figure 2a). From the top to bottom panels of
Figure 2c, we illustrate the QPO IMF and corresponding phase function,
instantaneous frequency and amplitude of the QPO component, respectively. As
one can see, the instantaneous frequency of QPO fluctuates between $\sim$1 and
$\sim$3 Hz over time.

Figure 4: The long-term combined power density spectra (PDS) of 30–100 (in
blue) and 100–200 (in red) keV energy bands. The shaded region highlights the
frequency range where the QPO fundamental is observed.

###  3.2 Recovery of QPO Waveforms above 120 keV

We perform the HHT analysis for each selected observations from HE/NaI in the
energy range of 30–120 keV. In such an energy band, HE shows a good
performance to provide high SNR light curves. These high SNR light curves are
suitable to be performed with the HHT analysis. Once the instantaneous phases,
$\phi(t)$, of the QPOs have been computed using the Hilbert transform as
plotted in Figure 2c, the QPO waveforms can be constructed by phase-folding
the light curves (also see Hu et al., 2014; Hsieh & Chou, 2020; Shui et al.,
2023a). We plot the reconstructed 30–120 keV QPO waveform in the top panel of
Figure 3. To enhance the SNR, waveforms from all selected observations have
been collectively combined, and the final combined waveform are normalized by
the mean value of the count rate. It is evident that the QPO waveform exhibits
an non-sinusoidal nature, which is characterized by a fast rise of $\sim 0.4$
cycles, followed by a gradual drop of $\sim 0.6$ cycles.

Utilizing the well-defined QPO phases derived from the HHT analysis, we extend
the analysis to perform phase-folding on high-energy light curves above 120
keV. It should be noted that, due to the low SNR in the high-energy band,
conducting the HHT analysis there is not feasible. Consequently, the phase-
folding of the high-energy light curves is carried out using the instantaneous
phase obtained from the HHT analysis in the low-energy band, assuming that the
high energy QPOs have the same frequency as that in the low energy band. From
the second to the bottom panels of Figure 3, the phase-folded light curves in
various energy bands above 120 keV are plotted. We use two methods to
determine the significance of the high energy QPO profiles, one is the
standard $\chi^{2}$-test method, another is the cross-correlation technique
described in Hou et al. (2022). It is evident that the QPO profile remains
distinctly identifiable in the energy bands between 120 and 200 keV. However,
signals dissipate in energy ranges exceeding 200 keV. The significance of the
QPO profile in the 150–170 keV energy band derived from the $\chi^{2}$-test
method is $5.2\sigma$, while the cross-correlation technique gives a
$8.1\sigma$ confidence level. As for the QPO profile in 170–200 keV, the
cross-correlation technique gives a $6.2\sigma$ confidence level (see Appendix
C for details). However, the $\chi^{2}$-test method gives only $2.1\sigma$ for
the 170–200 keV profile. This indicates that the cross-correlation method is
more advantageous since it takes the lower-energy QPO profile, which has a
very high SNR, as a template in searching for the higher-energy pulsed signals
(see also Hou et al., 2022). We also plot the PDS for high-energy light curves
(100–200 keV) to provide a comparison between the HHT phase-folding and
Fourier methods. Notably, no significant QPO peaks are observed in the long-
term combined PDS above 100 keV, as shown in Figure 4. In a simulation study
presented in Appendix B, we find that a small part of the noise component is
decomposed into the QPO IMF, resulting from mode mixing. However, the
simulation results show that if the QPO signals do not exist in the 170-200
keV energy band, the significance of the modulation from the noise is found to
be $<2\sigma$, which is inconsistent with the observational results ($\sim
6.2\sigma$). Conversely, when QPO signals are added to the high-energy light
curves, the significance of the modulation in the phase-folded light curve is
found to be $\sim 6-8\sigma$, consistent with the observational results. This
provides supports for the idea that the observed modulations of the phase-
folding profile in the energy ranges of 170–200 keV could indeed originate
from the QPO component, despite the presence of some mode mixing issues in the
VMD. (see Appendix B for details).

###  3.3 Properties of High Energy QPOs

It is crucial to highlight that the 170–200 keV QPO signal, as reported, marks
the highest-energy QPO detection for this source to date. To further study the
properties of the high-energy ($>$120 keV) QPO profiles, we analyze the phase-
folded waveforms using harmonically related cosine functions for fitting. From
Figure 3, it is evident that the QPO waveforms exhibit a non-sinusoidal
nature, particularly below 130 keV, which is consistent with the emergence of
the second harmonic in the PDS (see Figure 2b). We find that two cosine
functions provide a good description of the data. Furthermore, we subtract the
background contribution from the fractional amplitude, resulting in the final
fitting function of

$f(\phi)=1-\frac{S}{S+B}\\{A_{1}\cos{(\phi-\phi_{1})}+A_{2}\cos{\left[2(\phi-\phi_{2})\right]}\\},$
(7)

where $A_{1}$ and $A_{2}$ represent the amplitudes, $\phi_{1}$ and $\phi_{2}$
the phases of the cosine functions, and $S$ and $B$ denote the mean count
rates of the source and background, respectively. We perform phase-folding and
subsequent model fitting in 8 distinct energy bands of 30–40, 40–50, 50–70,
70–100, 100–130, 130–150, 150–170, and 170–200 keV. The analysis excludes
waveforms above 200 keV due to the insignificance of the QPO signal in these
bands. Notably, the parameters of the second harmonic component ($A_{2}$ and
$\phi_{2}$) are effectively constrained only at energies below 130 keV. This
suggests a notable attenuation of the second harmonic component above 130 keV,
prompting us to employ one cosine function for fitting in these energy bands.
Figure 5 presents the energy dependence of the best-fitting parameters. It is
clear that $\phi_{1}$ and $\phi_{2}$ remains relatively stable from 30 to 130
keV, whereas $\phi_{1}$ exhibits a decreasing trend with increasing energy
beyond 130 keV, indicating soft lags across different energy bands. Regarding
the amplitudes, while $A_{2}$ stays relatively constant, $A_{1}$ shows a
slight increase from 30 to 130 keV, followed by a decrease from $\sim 0.3$ to
$\sim 0.1$ above 130 keV.

Figure 5: The best-fitting phases (a) and amplitudes (b) of two harmonically
related cosine functions are plotted as a function of the photon energy.

##  4 Discussion

In this study, we have performed a comprehensive analysis of type-C QPOs in
the black hole binary MAXI J1535–571, utilizing data from the Insight-HXMT.
Based on the Hilbert-Huang transform (HHT) method, we have proposed a novel
method for detecting high-energy QPO signal. This allows us to successfully
reconstruct the QPO waveform above 100 keV through phase-folding of the light
curves, although no significant QPO signal has been discernible in the PDS
beyond this threshold. Notably, by combining QPO waveforms derived from seven
observations, we have achieved a highly confident detection of QPO profiles
above 170 keV at a significance level of $\sim 6.2\sigma$. This represents the
highest-energy detection of QPOs in this source reported to date. Furthermore,
this technique has the potential to extend the study of energy-dependent
behaviors of QPO amplitude and phase lag up to 200 keV.

The general method for detecting QPO signals from X-ray light curves is the
Fourier transform, which allows QPOs to be identified as a series of harmonic
peaks in the PDS. However, in the high-energy bands ($>$100 keV), the
limitations in count statistics present a challenge. Previous studies commonly
combined multiple observations to enhance the SNR of QPO peaks in the PDS (see
e.g. Ma et al., 2021; Yang et al., 2022). However, this approach assumes
similar QPO frequencies across the combined observations, which is not
typically the case during the HIMS of outburst rising phase of BHXRBs, where
the QPO frequency tends to increase rapidly (see e.g. Huang et al., 2018; Shui
et al., 2021, 2023b). If the QPO frequencies among combined observations show
significant differences, multiple fundamental QPO peaks would appear in the
long-term combined PDS in the low-energy band (see Figure 4). In such cases,
improving the SNR of higher energy QPO signals would be challenging, as the
QPO signals from different observations are distributed over a wide frequency
range. Although any method would face with the limitations in count statistics
in the high-energy bands, the HHT phase-folding method proposed in this work
effectively addresses the aforementioned issue encountered in the Fourier
method when combining different observations. In the phase-folding approach,
the well-defined QPO phases can be obtained in the low-energy band. This
indicates that, for any QPO frequency, a QPO cycle can be divided into given
phase bins between 0 and $2\pi$. Consequently, by stacking up counts within
the same phase bin from multiple observations, one can achieve a higher SNR
associated with the modulation of the count rate as function of the phase
(i.e. the QPO waveform). Both our observational and simulation studies
demonstrate that, provided the phase information of the QPO, the phase-folding
approach performs better in determining the high-energy QPO signal than the
Fourier method when the QPO frequency varies among different observations (see
Appendix B).

Previous studies by Huang et al. (2018) and Kong et al. (2020) have presented
the energy dependence of type-C QPO amplitude in MAXI J1535–571 across a wide
energy band ranging from 1 to 100 keV. They observed an increase in QPO RMS
amplitude with photon energy up to $\sim$20 keV, after which it plateaus
across all observations. Utilizing our HHT-based technique, we extend this
study to analyze QPO amplitude beyond 100 keV. As shown in Figure 5b, we note
a marginal increase in amplitude from 30 to 60 keV, followed by a plateau
between 60 and 120 keV, and a decreasing trend beyond 120 keV. The flat trend
of amplitude between 30 and 120 keV agrees with findings reported by Huang et
al. (2018) and Kong et al. (2020).

In broad-band spectral analyses of MAXI J1535–571, both Kong et al. (2020) and
Rodi et al. (2022), using Insight-HXMT and INTEGRAL/SPI data respectively,
reported findings indicative of a relatively low-temperature corona. Kong et
al. (2020) reported the cutoff energy of the power-law component within the
60–80 keV range, while Rodi et al. (2022) found the temperature of thermal
electrons in the corona to be $\sim 20$ keV during the HIMS. These
observations suggest that emissions above $\sim$100 keV are unlikely to be
corona-dominated. Furthermore, using INTEGRAL/SPI data, Rodi et al. (2022)
identified an additional high-energy component, which can be modeled by a
cutoff power-law and is observed to have a flux approximately six times over
the coronal flux in the 100–400 keV range (refer to Figure 6 in Rodi et al.
(2022)). Consequently, it is probable that the QPO amplitude above 100 keV is
predominantly influenced by this extra high-energy component. The different
evolutionary trends between the coronal and extra power-law fluxes, as
reported by Rodi et al. (2022), imply different origins for these two
components, with the high-energy component likely being from a jet. This is
further corroborated by the divergent energy-dependent properties of both
$\phi_{1}$ and $A_{1}$ in low- (30–130 keV) and high-energy (130–200 keV)
bands, as shown in Figure 5. Particularly, the emergence of soft lags above
130 keV agrees with predictions made by the jet precession model (Ma et al.,
2021, 2023; Shui et al., 2023a).

Using the phase function in 30–120 keV to fold the light curves above 120 keV
appears to introduce some systematic issues, as spectral analyses have
indicated different dominant components between the two energy bands. However,
our results demonstrate that the high-energy QPO waveform can be effectively
reconstructed using the 30–120 keV phase function, suggesting correlated QPO
phases across these energy bands. We propose that a jet-corona coupling
precession scenario could explain our findings, wherein the jet and corona are
coupled and precess with the same period, leading to correlated QPO signals
from the two locations. This picture is further supported by a recent general
relativistic magneto-hydrodynamic simulation (Liska et al., 2018).

Given that typical speeds at the base of the jet can already be relativistic
(McKinney, 2006), a fundamental concept in the jet precession model is that
emission along the jet axis is considerably stronger than off-axis emission
due to the Doppler boosting effect (see Ma et al., 2021). As a result, when
the jet precesses, the flux from it undergoes modulation. The velocity of the
jet material significantly influences the QPO variability. In particular,
higher jet speeds result in greater modulation amplitude, as there is a larger
difference in emission strength between on-axis and off-axis situations.
During the stage of the studied observations, it was reported to exist a
compact and continuous jet in MAXI J1535–571 from the radio band (Russell et
al., 2019). For a jet of this nature, the observed luminosity, $S_{\rm obs}$,
is described by

$S_{\rm obs}=S_{\rm int}\mathcal{D}^{\Gamma+2},$ (8)

where $S_{\rm obs}$ represents the observed luminosity, $S_{\rm int}$ is the
intrinsic luminosity of the jet, $\Gamma$ is the spectral index of the power-
law spectrum of the jet emission, and $\mathcal{D}$ is the relativistic
boosting factor given by

$\mathcal{D}=\frac{(1-\beta^{2})^{1/2}}{1-\beta\cos{\theta}},$ (9)

where $\beta=v/c$ is the jet speed in units of the speed of light, and
$\theta$ is the angle between the jet axis and the observer’s line of sight,
calculated as

$\displaystyle\cos{\theta}=$
$\displaystyle\sin{i}\cos{\Phi}\sin{\gamma}\cos{\omega}+\sin{i}\sin{\Phi}\sin{\gamma}\sin{\omega}$
$\displaystyle+\cos{\gamma}\cos{i},$ (10)

where $i$ is the inclination angle between the observer’s line of sight and
the black hole spin axis, $\Phi$ is the azimuth angle of the observer,
$\gamma$ is the angle between the jet axis and the black hole spin axis, and
$\omega$ is the precession phase angle, ranging from 0 to 2$\pi$ over the
precession period. (see Ma et al., 2021, for details). Setting the
inclination, $i$, to $57^{\circ}$ (Xu et al., 2018), and assuming $\gamma$ is
$5^{\circ}$ (see also Ma et al., 2021; Shui et al., 2023a) and $\Gamma=1.6$
which is the spectral index of the hard tail component (Rodi et al., 2022), we
employ the Markov Chain Monte Carlo (MCMC) method with a 50,000 sample length
to fit the QPO waveforms using the jet precession model. The model has two
free parameters: $\beta$ and $\Phi$. Given that the QPO radiative mechanism
above 130 keV is predominantly jet-driven, while below 130 keV it could be
affected by coronal emission, our fitting is only performed on the 130–150,
150–170, and 170–200 keV energy ranges. The fitting results are presented in
Figure 6. Based on the fitting with the jet precession model, we find
$\beta=0.44\pm 0.02$, $\beta=0.47\pm 0.06$ and $\beta=0.30\pm 0.07$ for
130–150, 150–170 and 170–200 keV energy bands, respectively. The value of
$\Phi$ depends on the reference point of the precession cycle, and the fitting
results of $\Phi=170^{\circ}.0\pm 3^{\circ}.7$, $\Phi=171^{\circ}.9\pm
8^{\circ}.9$ and $\Phi=143^{\circ}.4\pm 15^{\circ}.2$ for 130–150, 150–170 and
170–200 keV energy bands, respectively. The phase differences between these
energy bands represent soft lags. Based on the radio observations, Russell et
al. (2019) suggested an inclination of $\leqslant 45^{\circ}$. We also perform
fittings with the jet precession model, assuming $i=45^{\circ}$. The best-
fitting parameters for $\beta$ are $0.46\pm 0.02$, $0.49\pm 0.06$, and
$0.32\pm 0.07$ for the 130–150, 150–170, and 170–200 keV energy bands,
respectively, closely resembling the values obtained when $i=57^{\circ}$. The
uncertainties are reported at the $1\sigma$ confidence level.

Figure 6: The fitting results of the waveforms with the jet precession model
in the energy bands between 130 and 200 keV. Sequentially from left to right,
the fitting results for the energy bands of 130–150 keV, 150–170 keV, and
170–200 keV are depicted. All data points displayed have undergone background
subtraction. The red solid lines represent the jet precession model plotted
from the posterior distribution derived with the MCMC method.

##  5 Summary and Conclusion

A detailed analysis of type-C QPOs from the black hole binary MAXI J1535–571
with observations from Insight-HXMT is presented in this work. We propose a
novel method based on the Hilbert-Huang transform, which enables a robust
detection of type-C QPOs up to 200 keV, and consistently above 170 keV. This
method also has the potential to extend the study of the energy dependence of
QPO amplitude and phase lag up to 200 keV. Our findings reveal that the QPO
fundamental exhibits soft lags and a decrease in amplitude above $\sim 100$
keV. Given the reported detection of a hard tail in broad spectroscopy, we
proposed that the QPO properties above 100 keV are predominantly influenced by
the hard tail component, likely originating from a relativistic jet. Since the
high-energy QPO waveform can be well constructed using the phase function
derived from the low-energy light curves, a jet-corona coupling precession
scenario could explain our findings. Fitting the phase-folded profile with jet
precession model suggested a jet velocity of $\sim 0.4c$.

We are grateful to the anonymous referee for constructive comments that helped
us improve this paper. This research has made use of data obtained from the
High Energy Astrophysics Science Archive Research Center (HEASARC), provided
by NASA’s Goddard Space Flight Center, and the Insight-HXMT mission, a project
funded by China National Space Administration (CNSA) and the Chinese Academy
of Sciences (CAS). This work is supported by the National Key R&D Program of
China (2021YFA0718500) and the National Natural Science Foundation of China
under grants, U1838201, U1838202, 12173103, U2038101 and U1938103. This work
is partially supported by International Partnership Program of Chinese Academy
of Sciences (Grant No.113111KYSB20190020). Ling D. Kong is grateful for the
financial support provided bu the Sino-German (CSC-DAAD) Postdoc Scholarship
Program (57251553).

##  Appendix A VMD technique

VMD is a recently proposed technique designed to decompose a time series into
several intrinsic mode functions (Dragomiretskiy & Zosso, 2014). In the HHT
method, the original mode decomposition technique is the EMD (Huang et al.,
1998), which is an algorithm based on interpolation of local extrema. However,
this algorithm has a significant drawback: the problem of mode mixing, where
widely different scales could appear in a single IMF component, or a coherent
signal fragmented into separate parts that appear in more than one IMF
component. In comparison with the traditional EMD method, the VMD algorithm
theoretically reduces the impact of mode mixing by simultaneously decomposing
the input signal into a sum of IMFs while analytically determining the limited
center frequency and bandwidth.

In the VMD, the IMFs are represented as amplitude-modulated-frequency-
modulated (AM-FM) signals, given by $u_{k}(t)=A_{k}(t)\cos{(\phi_{k}(t))}$,
where the phase for the $k$th mode, $\phi_{k}(t)$, is a non-decreasing
function, and the envelope $A_{k}(t)$ is non-negative ($A_{k}(t)\geq 0$). When
applied to an input signal $f(t)$, the VMD initializes the set of all modes as
$\\{u_{k}\\}=\\{u_{1},...,u_{K}\\}$ along with their corresponding center
frequencies $\\{\omega_{k}\\}=\\{\omega_{1},...,\omega_{K}\\}$. The goal is to
minimize the sum of bandwidths of IMFs by solving the optimization problem
(1). To estimate the bandwidth of each mode (BKk), the VMD employs several
signal processing tools. Details on these signal processing tools can be found
in Dragomiretskiy & Zosso (2014). The bandwidth of each mode is estimated
using $H^{1}$ Gaussian smoothness of the shifted signal, i.e. the $L^{2}$ norm
of the gradient:

${\rm BW}_{k}=\left|\left|\partial_{t}\left[\left(\delta(t)+\frac{j}{\pi
t}\right)*u_{k}(t)\right]e^{-j\omega_{k}t}\right|\right|_{2},$ (A1)

where $\delta$ is the Dirac function, and ‘$*$’ denotes convolution. The term
enclosed in the square brackets represents the analytic signal of $u_{k}(t)$.
With Lagrangian multipliers, $\lambda$, the augmented Lagrangian,
$\mathcal{L}$, can be expressed as

$\displaystyle\mathcal{L}(\\{u_{k}\\},\\{\omega_{k}\\},\lambda)=\alpha\sum_{k=1}^{K}{\rm
BW}_{k}^{2}+\left|\left|f(t)-\sum_{k=1}^{K}u_{k}(t)\right|\right|_{2}^{2}+\left\langle\lambda(t),f(t)-\sum_{k=1}^{K}u_{k}(t)\right\rangle.$
(A2)

The solution to the original minimization problem (1) is now achieved as the
saddle point of the augmented Lagrangian, which can be found using the
alternate direction multiplier method (ADMM). Then all the modes gained from
solutions in Fourier domain are written as

$\hat{u}_{k}(\omega)=\frac{\hat{f}(\omega)-\sum_{i\neq
k}{\hat{u_{i}}(\omega)}+\frac{\hat{\lambda}(\omega)}{2}}{1+2\alpha(\omega-\omega_{k})^{2}},$
(A3)

where $\hat{f}(\omega)$ and $\hat{u}_{k}(\omega)$ represent the Fourier
transform of the signal $f(t)$ and mode $u_{k}(t)$, respectively, $\omega$ is
the Fourier frequency. The VMD algorithm iteratively refines the modes
$\\{u_{k}\\}$ and their corresponding center frequencies $\\{\omega_{k}\\}$
until they stabilize. In this algorithm, the center frequency $\omega_{k}$ is
determined within the Fourier domain by the formula:

$\omega_{k}=\frac{\int_{0}^{\infty}\omega\left|\hat{u}_{k}(\omega)\right|^{2}d\omega}{\int_{0}^{\infty}\left|\hat{u}_{k}(\omega)\right|^{2}d\omega}.$
(A4)

This approach effectively extends Wiener filtering to accommodate adaptive
filtering across multiple frequency bands. For a comprehensive exploration of
this technique and the constrained variational optimization problem, see
Dragomiretskiy & Zosso (2014).

##  Appendix B Simulation Analysis

Figure 7: Power density spectrum of the HXMT/HE observation P011453500201 in
the 30–120 keV energy band fitted with a model consists of a broken power law
function (plotted in red) and two Lorentzian functions (plotted in blue). In
the plotted power density spectrum, Poisson noise has been subtracted.

Figure 8: Left: simulated light curve composed of a red noise and a QPO
component. In this simulated signal, the Poisson noise is also considered.
Right: Fourier power spectra that are produced from the original total light
curve (in black), simulated QPO component (in grey), and decomposed QPO
component (in blue), respectively.

Figure 9: The top panels respectively compare the time variations of the QPO
count rate and instantaneous phase between the decomposed QPO IMF (in blue)
and synthetic QPO signal (in grey) in a representative 20-second time range.
The two middle panels present the parameter correlations between decomposed
and synthetic QPO signals using all time bins. The two bottom panels show the
histograms of count rate difference and phase difference between the simulated
QPO and decomposed QPO signals from all time bins. The orange dashed lines
displayed in the middle and bottom panels represent the case that the two
signals are in perfect alignment.

To assess the reliability of the QPO phase-folding analysis results, we
conduct a comprehensive simulation analysis for our HHT-based method to test
its robustness. This simulation analysis consists of three parts: (1) testing
the performance of the VMD algorithm to verify that a simulated QPO, as strong
as the observational one, is indeed decomposed accurately; (2) exploring how
the VMD algorithm deals with a “pure noise” signal; and (3) addressing a
potential issue that may arise from (1) and (2) in the phase-folding analysis.

While it has now become standard practice to fit power spectra in X-ray
binaries with a combination of both broad and peaked components modeled as
Lorentzians (see e.g. Belloni et al., 2002; Yang et al., 2022; Zhou et al.,
2022; Shui et al., 2023b), there is currently no physical backing behind the
use of Lorentzian components. Additionally, we find the continuum of the power
spectrum from obsID P011453500201 in 30–120 keV can be well fitted with a
broken power law function (see Figure 7). Consequently, we individually model
the red noise with a broken power law function and QPO harmonics with the
Lorentzian functions. Using the function shapes derived from the modeling of
the observed power spectrum, we apply the algorithm from Timmer & Koenig
(1995) to generate 1684-s light curves with a 0.05-s time resolution for both
red noise and QPO components. In this algorithm, the power spectral amplitude
at each frequency is randomly drawn from a $\chi^{2}$ distribution, and not
fixed at the amplitude of the underlying power spectrum. Additionally, the
Fourier phases given to the signal by this algorithm is random. For the broken
power law funtion of the red noise spectrum, we set the two power law indices
at 0.28 and 1.59, respectively, the break frequency at 5.61 Hz, and the
fractional rms amplitude at 30%. As for the Lorentzian function of the QPO
component, we set the centroid frequency at 1.78 Hz, the FWHM at 0.12 Hz, and
the fractional rms amplitude at 15%, respectively. The simulations of the
light curves are performed using
Stingray555https://github.com/StingraySoftware/stingray/tree/v1.1.2, which is
an open-source Python package (Huppenkothen et al., 2019). In this study, we
focus on the QPO fundamental and therefore do not include the harmonic
component in the simulation analysis. The total light curve with an average
count rate of 600 cts s-1 consists of the red noise, QPO and Poisson noise,
which is presented in the left panel of Figure 8. Subsequently, the VMD
technique is used to decompose this light curve. It is also important to note
that Poisson noise is the observational errors caused by statistical
fluctuations in photon counting, while red noise is related to the intrinsic
stochastic behavior of the source. Therefore, Poisson noise is specifically
incorporated by Poisson sampling of the variable source flux which includes
both the red noise and QPO components. The right panel of Figure 8 shows the
power spectra of the simulated total (in black), QPO (in grey) light curves
and decomposed QPO IMF (in blue). As one can see, the power spectrum of the
decomposed QPO IMF is well consistent with that of the simulated QPO signal.
This demonstrates that the VMD can precisely recover the essential features of
the QPO in the Fourier domain, including the centroid frequency, width and
amplitude. We then perform additional Hilbert transform analyses on both the
decomposed and simulated QPO light curves. This is followed by comparing the
trends in QPO count rate and instantaneous phase between the two sets of light
curves. The top panels of Figure 9 present a comparison of the time variations
in QPO count rate and instantaneous phases between the decomposed (in blue)
and simulated (in grey) QPO signals in a representative 20-s time range. The
middle panels present the parameter correlations between the two signals,
incorporating all time bins. In addition, the two bottom panels show the
histograms of count rate difference and phase difference between the simulated
QPO and decomposed QPO signals, respectively. The orange dashed lines
displayed in the middle and bottom panels represent the case that the two
signals are in perfect alignment. We find the structures of the decomposed QPO
are highly consistent with those of the simulated QPO signal, but with local
fluctuations, possibly due to mode mixing between the QPO and noise
components, or the statistical fluctuations. Specifically, in the right middle
panel of Figure 9, a small number of data points are distributed in the upper
left and lower right of the plot, corresponding to the phase difference
distributed around $-2\pi$ and $2\pi$ in the right bottom panel. Since the
phase of the QPO varies from $-\pi$ to $\pi$, this phenomenon reflects that
the phases of the simulated and decomposed QPO signals in a small part of time
bins have local shifts.

Figure 10: Power density spectra of the simulated pure Poisson noise (in
black) and the corresponding VMD IMFs (in colors). The decomposed IMFs are
evenly distributed in frequencies, suggesting that VMD is essentially a set of
adaptive Wiener filter bank capable of separating modes with different center
frequencies.

The VMD technique is also applied to a pure Poisson noise signal. In Figure
10, we plot the power spectra of the Poisson noise and corresponding VMD IMFs.
The result shows that the decomposed IMFs are evenly distributed in
frequencies, suggesting that VMD is essentially a set of adaptive Wiener
filter bank capable of separating modes with different center frequencies. The
equivalent filter bank properties of the VMD have been extensively studied in
literature (see e.g. Wang et al., 2015). From another perspective, the VMD
algorithm is capable of generating the output IMFs for any given signal, even
if the signal is pure noise. This indicates that the IMF corresponding to the
QPO component might be influenced by the noises, potentially leading to local
variations between the real and decomposed QPO signals observed in Figure 9.

Figure 11: Phase-folding profiles of the simulated total (in blue) and noise
(in red) light curves. The phase-folding is performed by using the
instantaneous phase function of the decomposed QPO IMFs and by combining the
seven simulated light curves (see the corresponding text for details). Figure
12: Power density spectra of the simulated light curves in the energy ranges
of 150–170 (in orange) and 170–200 (in blue) keV.

To explore the potential impact of the local variations between simulated and
decomposed QPO light curves presented in Figure 9 on the phase-folded light
curve, we conduct the phase-folding on the initially simulated total and QPO-
subtracted light curves, utilizing the instantaneous phase of the decomposed
QPO IMF. Here, the QPO-subtracted light curve refers to the red noise with its
Poisson noise. In our observational study, we combined data from seven
different observations, each with varying QPO frequencies and exposures. Our
simulation analysis is designed to closely replicate the observational
analysis, taking into account these variations in QPO frequencies and
exposures as accurately as possible. Therefore, we simulate seven sets of
light curves in which the QPO frequency and exposure times are set to be the
same as those in the corresponding observations. For each total light curve,
we apply the VMD method and the Hilbert transform. Subsequently, the phase-
folded total and QPO-subtracted light curves are obtained by phase-folding
across the seven light curves using the QPO instantaneous phases obtained from
the HHT analysis. As shown in Figure 11, the phase-folding profile of the
total light curve exhibits a significant modulation, which is defined as the
QPO waveform. The amplitude of the modulation is estimated to be $\sim 0.21$,
with a corresponding rms of $\sim 14.9$%, consistent with the initially set
value in the simulation (15%). However, a relatively weaker modulation (with
an amplitude of $\sim 0.04$) of phase-folded QPO-subtracted light curve is
also observed. This indicates that a small portion of the red noise is
decomposed into the QPO IMF, resulting from the issue of mode mixing.

As the phase-folded noise light curve simulated in the low-energy band can
exhibit modulations due to mode mixing (see Figure 11), it is important to
confirm that the observed modulations in the high-energy bands (see Figure 3)
are indeed from the QPO component and not from the red noise. To achieve this,
we simulate high-energy light curves for the seven observations based on the
conditions of the observational data. Specifically, in the 150–170 and 170–200
keV energy bands, the average total count rates of the observed light curves
are $\sim 32$ and $47$ cts s-1, respectively, with the net count rates of
$\sim 2$ cts s-1 in both energy bands. For simplicity, we assume that the
light curves in different energy bands for both the QPO and red noise are
fully coherent. This indicates that the light-curves of both the QPO and red
noise have consistent phase functions across different energy bands.
Additionally, we assume that the fractional rms of the red noise and the QPO
remains constant with photon energy. Therefore, to simulate red noise and QPO
light curves in the two high-energy bands, we scale the absolute amplitude of
the original simulated red noise and QPO light curves in the 30–120 keV energy
band by a factor of 1/300 ($\sim 2$ cts s-1/$\sim 600$ cts s-1). Additionally,
in these simulated total light curves, Poisson noise is also incorporated.

Following the analysis presented in Section 3, for each energy band, we
combine these seven light curves to generate the PDS. Similar to the
observational findings presented in Figure 4, the combined power spectra of
the simulated light curves in the high-energy bands are dominant with Poisson
noise, resulting in QPO peaks that are too weak to be robustly detected (see
Figure 12). To rigorously evaluate the phase-folding profiles under the impact
of Poisson noise, we independently simulate $10^{4}$ light curves for each of
the two energy bands. Then, the phase-folding is performed on the simulated
total and QPO-subtracted light curves for both energy bands. It should be
noted that, due to the low SNR in the high-energy band, conducting the HHT
analysis there is not feasible. Consequently, the phase-folding of the high-
energy light curves is carried out using the instantaneous phase obtained from
the HHT analysis in the low-energy band.

The upper panels of Figure 13 present phase-folding profiles using 1000
samples from the simulated light curves in the two energy bands. It is evident
that in both energy bands, the phase-folding profiles of the total light
curves exhibit significant modulations, while those of the QPO-subtracted
light curves are relatively weaker, likely being dominated by the random
fluctuations of Poisson noise. Based on the cross-correlation technique
described in Appendix C, for each sample, we calculate the maximum cross-
correlations between phase-folding profiles of the total, QPO-subtracted, pure
Poisson noise light curves in the high-energy bands and the reference profile.
The reference profile is the phase-folded light curve in the low-energy band,
plotted in blue in Figure 11.

As shown in the bottom panels of Figure 13, the distributions of maximum
cross-correlations of the pure Poisson noise and QPO-subtracted signals
exhibit significant overlap in both energy bands. This demonstrates the QPO-
subtracted signals in the two energy bands are dominated by the Poisson noise.
Therefore, the modulations in the phase-folding profiles of the QPO-subtracted
signals in the high-energy bands could only be determined at the confidence
level of $<2\sigma$, which is inconsistent with our observational results
($\sim$6–8 $\sigma$). Moreover, in a real case, the red noise light curves
from two energy bands are not perfectly correlated, indicating that the
modulations of the QPO-subtracted signals in the high energy bands could be
even weaker than the simulated ones. However, the distributions of the total
signals (with the QPO component) are significantly different from those of the
pure Poisson and QPO-subtracted signals, indicating that the modulations of
total light curves in the high energy bands could be detected at the $\sim$5–8
$\sigma$ confidence level, consistent with the observational results.

In the simulation analysis above, we assume that the fractional rms of QPO and
red noise does not vary with photon energy. However, in reality, the
fractional rms of both signals may indeed change with photon energy.
Therefore, we conduct tests to examine how the amplitude of the QPO and red
noise affects the phase-folding method performed in the 150–170 keV and
170–200 keV energy bands. Firstly, we keep the fractional rms of the red noise
constant at the value observed in the 30–120 keV energy band, while varying
the QPO fractional rms. The top two panels of Figure 14 present the dependence
of the modulation significance on the QPO rms. The significance is determined
with the Monte Carlo method described in Appendix C. It is evident that the
modulation significance of the total profiles increases with the QPO rms,
whereas that of the QPO-subtracted profiles remain relatively constant
(0–4$\sigma$). When the QPO fractional rms is lower than $\sim 8\%$, the
modulation of the total profiles cannot be robustly determined ($<5\sigma$).
Subsequently, we fix the fractional rms of the QPO signal at the value of that
in the 30–120 keV energy band, and vary the fractional rms of the red noise.
As shown in the two bottom panels of Figure 14, both the total and QPO-
subtracted profiles exhibit modulations that can be determined with higher
confidence levels when the fractional rms of the red noise is higher. However,
for the cases that the fractional rms of the red noise is lower than 60%, the
modulations of the QPO-subtracted profiles could only be determined at the
confidence level of $<5\sigma$, which is inconsistent with our observational
results ($\sim$6–8 $\sigma$). Moreover, previous studies on the red noise in
BHXRBs have indicated that its fractional rms generally decreases with photon
energy (see e.g. Bu et al., 2021; Yang et al., 2022). This suggests that the
fractional rms of the red noise is likely to be lower than 30% in the 150–170
keV and 170–200 keV energy bands, and could not significantly influence the
determination of the intrinsic QPO modulation.

In our simulation analysis, we model the red noise component using a broken
power law function. We have also conducted the simulation analysis using the
zero-frequency Lorentzian function to model the red noise, and obtained
similar results. Our simulation analyses provide supports for the idea that
the observed modulations of the phase-folding profile in the energy ranges of
150–170 and 170–200 keV could indeed originate from the QPO component, despite
the presence of some mode mixing issues in the VMD.

Figure 13: The top two panels respectively present phase-folding profiles of
total (in blue) and QPO-subtracted (in red) light curves sampled from 104
simulated light curves in the energy ranges of 150–170 and 170–200 keV. The
bottom two panels present the distributions of the maximum cross-correlations
between phase-folding profiles of the simulated total (in blue), QPO-
subtracted (in red), pure Poisson noise (in grey) light curves in the two high
energy bands and the profile of the total light curve in the low energy band
presented in Figure 11.

Figure 14: The top two panels respectively present relationships between QPO
fractional rms and significance of modulations in the total (in blue) and QPO-
subtracted (in red) profiles in the energy ranges of 150–170 and 170–200 keV,
where the red noise fractional rms is fixed at 30%. The bottom two panels
present relationships between red noise fractional rms and significance of
modulations in the total and QPO-subtracted profiles, where the QPO fractional
rms is fixed at 15%. The purple dotted lines plotted in four panels
represented that the modulation in profiles can be determined at a $5\sigma$
confidence level.

##  Appendix C Determine the high-energy QPO waveforms

In this section, we present the determination of the high-energy QPO waveforms
using a cross-correlation technique (see Hou et al., 2022, for more details).
In a brief, this method first calculates the cross correlation of the QPO
waveform in a given energy band with that in the 30–120 keV energy range, and
then uses a Monte Carlo method to simulate the cross-correlation distribution
between two noncorrelated Poisson sampled profiles. Appendix Figure 15 shows
the cross correlation of QPO waveforms and the distribution of the cross
correlation from simulated waweforms for energy ranges of 120–130 and 170–200
keV as examples. The number of bins per period is chosen as 60, but we have
tested different binning of 30, 40, 50, 70, all resulting in similar
significance of the QPO waveform in a given energy band. Using this technique,
we detect the high-energy QPO waveform above 170 keV up to 200 keV with a
significance of 6.2$\sigma$.

Figure 15: Examples of cross correlation of 120–130 and 170–200 keV QPO
waveforms with the 30–120 keV QPO waveform. The upper panels plot the cross
correlation vs. phase lag for a specified waveform with the 30–120 keV QPO
waveform. The lower panels show the cross-correlation distribution from
$10^{4}$ Monte Carlo simulations of two non-correlated waveforms. The dashed
vertical lines plotted in the lower panels indicate the peak cross-correlation
values as observed in the upper panels.

## References

  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 813, L41, doi: 10.1088/2041-8205/813/2/L41
  * Axelsson & Done (2016) Axelsson, M., & Done, C. 2016, MNRAS, 458, 1778, doi: 10.1093/mnras/stw464
  * Axelsson et al. (2014) Axelsson, M., Done, C., & Hjalmarsdotter, L. 2014, MNRAS, 438, 657, doi: 10.1093/mnras/stt2236
  * Baglio et al. (2018) Baglio, M. C., Russell, D. M., Casella, P., et al. 2018, ApJ, 867, 114, doi: 10.3847/1538-4357/aae532
  * Banerjee et al. (2023) Banerjee, A., Sharma, A., Mandal, A., et al. 2023, MNRAS, 523, L52, doi: 10.1093/mnrasl/slad057
  * Bellavita et al. (2022) Bellavita, C., García, F., Méndez, M., & Karpouzas, K. 2022, MNRAS, 515, 2099, doi: 10.1093/mnras/stac1922
  * Belloni et al. (2005) Belloni, T., Homan, J., Casella, P., et al. 2005, A&A, 440, 207, doi: 10.1051/0004-6361:20042457
  * Belloni et al. (2002) Belloni, T., Psaltis, D., & van der Klis, M. 2002, ApJ, 572, 392, doi: 10.1086/340290
  * Belloni & Stella (2014) Belloni, T. M., & Stella, L. 2014, Space Sci. Rev., 183, 43, doi: 10.1007/s11214-014-0076-0
  * Bhargava et al. (2019) Bhargava, Y., Belloni, T., Bhattacharya, D., & Misra, R. 2019, MNRAS, 488, 720, doi: 10.1093/mnras/stz1774
  * Bu et al. (2021) Bu, Q. C., Zhang, S. N., Santangelo, A., et al. 2021, ApJ, 919, 92, doi: 10.3847/1538-4357/ac11f5
  * Camp et al. (2007) Camp, J. B., Cannizzo, J. K., & Numata, K. 2007, Phys. Rev. D, 75, 061101, doi: 10.1103/PhysRevD.75.061101
  * Cao et al. (2020) Cao, X., Jiang, W., Meng, B., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249504, doi: 10.1007/s11433-019-1506-1
  * Carvalho et al. (2020) Carvalho, V. R., Moraes, M. F., Braga, A. P., & Mendes, E. M. 2020, Biomedical Signal Processing and Control, 62, 102073, doi: https://doi.org/10.1016/j.bspc.2020.102073
  * Casella et al. (2005) Casella, P., Belloni, T., & Stella, L. 2005, ApJ, 629, 403, doi: 10.1086/431174
  * Chauhan et al. (2019) Chauhan, J., Miller-Jones, J. C. A., Anderson, G. E., et al. 2019, MNRAS, 488, L129, doi: 10.1093/mnrasl/slz113
  * Chen et al. (2020) Chen, Y., Cui, W., Li, W., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249505, doi: 10.1007/s11433-019-1469-5
  * de Ruiter et al. (2019) de Ruiter, I., van den Eijnden, J., Ingram, A., & Uttley, P. 2019, MNRAS, 485, 3834, doi: 10.1093/mnras/stz665
  * Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&A Rev., 15, 1, doi: 10.1007/s00159-007-0006-1
  * Dragomiretskiy & Zosso (2014) Dragomiretskiy, K., & Zosso, D. 2014, IEEE Transactions on Signal Processing, 62, 531, doi: 10.1109/TSP.2013.2288675
  * Guo et al. (2020) Guo, C.-C., Liao, J.-Y., Zhang, S., et al. 2020, Journal of High Energy Astrophysics, 27, 44, doi: 10.1016/j.jheap.2020.02.008
  * Homan & Belloni (2005) Homan, J., & Belloni, T. 2005, Ap&SS, 300, 107, doi: 10.1007/s10509-005-1197-4
  * Homan et al. (2020) Homan, J., Bright, J., Motta, S. E., et al. 2020, ApJ, 891, L29, doi: 10.3847/2041-8213/ab7932
  * Hou et al. (2022) Hou, X., Ge, M. Y., Ji, L., et al. 2022, ApJ, 938, 149, doi: 10.3847/1538-4357/ac8c93
  * Hsieh & Chou (2020) Hsieh, H.-E., & Chou, Y. 2020, ApJ, 900, 116, doi: 10.3847/1538-4357/abacbd
  * Hu et al. (2011) Hu, C.-P., Chou, Y., Wu, M.-C., Yang, T.-C., & Su, Y.-H. 2011, ApJ, 740, 67, doi: 10.1088/0004-637X/740/2/67
  * Hu et al. (2014) Hu, C.-P., Chou, Y., Yang, T.-C., & Su, Y.-H. 2014, ApJ, 788, 31, doi: 10.1088/0004-637X/788/1/31
  * Huang & Wu (2008) Huang, N. E., & Wu, Z. 2008, Reviews of Geophysics, 46, RG2006, doi: 10.1029/2007RG000228
  * Huang et al. (1998) Huang, N. E., Shen, Z., Long, S. R., et al. 1998, Proceedings of the Royal Society of London Series A, 454, 903, doi: 10.1098/rspa.1998.0193
  * Huang et al. (2018) Huang, Y., Qu, J. L., Zhang, S. N., et al. 2018, ApJ, 866, 122, doi: 10.3847/1538-4357/aade4c
  * Huppenkothen et al. (2019) Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, ApJ, 881, 39, doi: 10.3847/1538-4357/ab258d
  * Ingram & Done (2011) Ingram, A., & Done, C. 2011, MNRAS, 415, 2323, doi: 10.1111/j.1365-2966.2011.18860.x
  * Ingram et al. (2009) Ingram, A., Done, C., & Fragile, P. C. 2009, MNRAS, 397, L101, doi: 10.1111/j.1745-3933.2009.00693.x
  * Ingram & van der Klis (2015) Ingram, A., & van der Klis, M. 2015, MNRAS, 446, 3516, doi: 10.1093/mnras/stu2373
  * Ingram et al. (2017) Ingram, A., van der Klis, M., Middleton, M., Altamirano, D., & Uttley, P. 2017, MNRAS, 464, 2979, doi: 10.1093/mnras/stw2581
  * Ingram et al. (2016) Ingram, A., van der Klis, M., Middleton, M., et al. 2016, MNRAS, 461, 1967, doi: 10.1093/mnras/stw1245
  * Ingram & Motta (2019) Ingram, A. R., & Motta, S. E. 2019, New A Rev., 85, 101524, doi: 10.1016/j.newar.2020.101524
  * Kato (1990) Kato, S. 1990, PASJ, 42, 99
  * Kennea et al. (2017) Kennea, J. A., Evans, P. A., Beardmore, A. P., et al. 2017, The Astronomer’s Telegram, 10700, 1
  * Kong et al. (2020) Kong, L. D., Zhang, S., Chen, Y. P., et al. 2020, Journal of High Energy Astrophysics, 25, 29, doi: 10.1016/j.jheap.2020.01.003
  * Lachowicz & Done (2010) Lachowicz, P., & Done, C. 2010, A&A, 515, A65, doi: 10.1051/0004-6361/200913144
  * Liao et al. (2020a) Liao, J.-Y., Zhang, S., Lu, X.-F., et al. 2020a, Journal of High Energy Astrophysics, 27, 14, doi: 10.1016/j.jheap.2020.04.002
  * Liao et al. (2020b) Liao, J.-Y., Zhang, S., Chen, Y., et al. 2020b, Journal of High Energy Astrophysics, 27, 24, doi: 10.1016/j.jheap.2020.02.010
  * Liska et al. (2018) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81, doi: 10.1093/mnrasl/slx174
  * Liu et al. (2020) Liu, C., Zhang, Y., Li, X., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249503, doi: 10.1007/s11433-019-1486-x
  * Liu et al. (2021) Liu, H.-X., Huang, Y., Xiao, G.-C., et al. 2021, Research in Astronomy and Astrophysics, 21, 070, doi: 10.1088/1674-4527/21/3/70
  * Ma et al. (2021) Ma, X., Tao, L., Zhang, S.-N., et al. 2021, Nature Astronomy, 5, 94, doi: 10.1038/s41550-020-1192-2
  * Ma et al. (2023) Ma, X., Zhang, L., Tao, L., et al. 2023, ApJ, 948, 116, doi: 10.3847/1538-4357/acc4c3
  * McKinney (2006) McKinney, J. C. 2006, MNRAS, 368, 1561, doi: 10.1111/j.1365-2966.2006.10256.x
  * Méndez et al. (2022) Méndez, M., Karpouzas, K., García, F., et al. 2022, Nature Astronomy, 6, 577, doi: 10.1038/s41550-022-01617-y
  * Mereminskiy et al. (2018) Mereminskiy, I. A., Grebenev, S. A., Prosvetov, A. V., & Semena, A. N. 2018, Astronomy Letters, 44, 378, doi: 10.1134/S106377371806004X
  * Mineshige et al. (1994) Mineshige, S., Takeuchi, M., & Nishimori, H. 1994, ApJ, 435, L125, doi: 10.1086/187610
  * Molteni et al. (1996) Molteni, D., Sponholz, H., & Chakrabarti, S. K. 1996, ApJ, 457, 805, doi: 10.1086/176775
  * Morgan et al. (1997) Morgan, E. H., Remillard, R. A., & Greiner, J. 1997, ApJ, 482, 993, doi: 10.1086/304191
  * Motta et al. (2015) Motta, S. E., Casella, P., Henze, M., et al. 2015, MNRAS, 447, 2059, doi: 10.1093/mnras/stu2579
  * Nakahira et al. (2018) Nakahira, S., Shidatsu, M., Makishima, K., et al. 2018, PASJ, 70, 95, doi: 10.1093/pasj/psy093
  * Nathan et al. (2022) Nathan, E., Ingram, A., Homan, J., et al. 2022, MNRAS, 511, 255, doi: 10.1093/mnras/stab3803
  * Nazari & Sakhaei (2018) Nazari, M., & Sakhaei, S. M. 2018, IEEE Journal of Biomedical and Health Informatics, 22, 1059, doi: 10.1109/JBHI.2017.2734074
  * Negoro et al. (2017) Negoro, H., Ishikawa, M., Ueno, S., et al. 2017, The Astronomer’s Telegram, 10699, 1
  * Nolan et al. (1981) Nolan, P. L., Gruber, D. E., Matteson, J. L., et al. 1981, ApJ, 246, 494, doi: 10.1086/158949
  * Parikh et al. (2019) Parikh, A. S., Russell, T. D., Wijnands, R., et al. 2019, ApJ, 878, L28, doi: 10.3847/2041-8213/ab2636
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Rodi et al. (2022) Rodi, J., Jourdain, E., & Roques, J. P. 2022, ApJ, 935, 25, doi: 10.3847/1538-4357/ac7fff
  * Russell et al. (2019) Russell, T. D., Tetarenko, A. J., Miller-Jones, J. C. A., et al. 2019, ApJ, 883, 198, doi: 10.3847/1538-4357/ab3d36
  * Sandrinelli et al. (2014) Sandrinelli, A., Covino, S., & Treves, A. 2014, ApJ, 793, L1, doi: 10.1088/2041-8205/793/1/L1
  * Scargle et al. (1993) Scargle, J. D., Steiman-Cameron, T., Young, K., et al. 1993, ApJ, 411, L91, doi: 10.1086/186920
  * Shui et al. (2021) Shui, Q. C., Yin, H. X., Zhang, S., et al. 2021, MNRAS, 508, 287, doi: 10.1093/mnras/stab2521
  * Shui et al. (2023a) Shui, Q. C., Zhang, S., Zhang, S. N., et al. 2023a, ApJ, 957, 84, doi: 10.3847/1538-4357/acfc42
  * Shui et al. (2023b) Shui, Q. C., Zhang, S., Chen, Y. P., et al. 2023b, ApJ, 943, 165, doi: 10.3847/1538-4357/aca7b8
  * Sobolewska & Życki (2006) Sobolewska, M. A., & Życki, P. T. 2006, MNRAS, 370, 405, doi: 10.1111/j.1365-2966.2006.10489.x
  * Stella et al. (1999) Stella, L., Vietri, M., & Morsink, S. M. 1999, ApJ, 524, L63, doi: 10.1086/312291
  * Stevens & Uttley (2016) Stevens, A. L., & Uttley, P. 2016, MNRAS, 460, 2796, doi: 10.1093/mnras/stw1093
  * Stevens et al. (2018) Stevens, A. L., Uttley, P., Altamirano, D., et al. 2018, ApJ, 865, L15, doi: 10.3847/2041-8213/aae1a4
  * Stiele & Kong (2018) Stiele, H., & Kong, A. K. H. 2018, ApJ, 868, 71, doi: 10.3847/1538-4357/aae7d3
  * Su et al. (2015) Su, Y.-H., Chou, Y., Hu, C.-P., & Yang, T.-C. 2015, ApJ, 815, 74, doi: 10.1088/0004-637X/815/1/74
  * Tagger & Pellat (1999) Tagger, M., & Pellat, R. 1999, A&A, 349, 1003
  * Tao et al. (2018) Tao, L., Chen, Y., Güngör, C., et al. 2018, MNRAS, 480, 4443, doi: 10.1093/mnras/sty2157
  * Terrell (1972) Terrell, N. James, J. 1972, ApJ, 174, L35, doi: 10.1086/180944
  * Tetarenko et al. (2017) Tetarenko, A. J., Russell, T. D., Miller-Jones, J. C. A., Sivakoff, G. R., & Jacpot Xrb Collaboration. 2017, The Astronomer’s Telegram, 10745, 1
  * Timmer & Koenig (1995) Timmer, J., & Koenig, M. 1995, A&A, 300, 707
  * van den Eijnden et al. (2017) van den Eijnden, J., Ingram, A., Uttley, P., et al. 2017, MNRAS, 464, 2643, doi: 10.1093/mnras/stw2634
  * van der Klis (1989) van der Klis, M. 1989, ARA&A, 27, 517, doi: 10.1146/annurev.aa.27.090189.002505
  * Wagoner (1999) Wagoner, R. V. 1999, Phys. Rep., 311, 259, doi: 10.1016/S0370-1573(98)00104-5
  * Wang et al. (2015) Wang, Y., Markert, R., Xiang, J., & Zheng, W. 2015, Mechanical Systems and Signal Processing, 60-61, 243, doi: https://doi.org/10.1016/j.ymssp.2015.02.020
  * Wijnands et al. (1999) Wijnands, R., Homan, J., & van der Klis, M. 1999, ApJ, 526, L33, doi: 10.1086/312365
  * Xu et al. (2018) Xu, Y., Harrison, F. A., García, J. A., et al. 2018, ApJ, 852, L34, doi: 10.3847/2041-8213/aaa4b2
  * Yang et al. (2022) Yang, Z.-X., Zhang, L., Bu, Q.-C., et al. 2022, ApJ, 932, 7, doi: 10.3847/1538-4357/ac63af
  * Yu et al. (2023a) Yu, W., Bu, Q.-C., Yang, Z.-X., et al. 2023a, ApJ, 951, 130, doi: 10.3847/1538-4357/acd756
  * Yu et al. (2023b) Yu, W., Bu, Q.-C., Liu, H.-X., et al. 2023b, ApJ, 953, 191, doi: 10.3847/1538-4357/acd9a2
  * Zhang et al. (2017) Zhang, L., Wang, Y., Méndez, M., et al. 2017, ApJ, 845, 143, doi: 10.3847/1538-4357/aa8138
  * Zhang et al. (2020a) Zhang, L., Méndez, M., Altamirano, D., et al. 2020a, MNRAS, 494, 1375, doi: 10.1093/mnras/staa797
  * Zhang et al. (2023) Zhang, P., Soria, R., Zhang, S., et al. 2023, A&A, 677, A178, doi: 10.1051/0004-6361/202346309
  * Zhang et al. (2014) Zhang, S., Lu, F. J., Zhang, S. N., & Li, T. P. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9144, Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray, ed. T. Takahashi, J.-W. A. den Herder, & M. Bautz, 914421, doi: 10.1117/12.2054144
  * Zhang et al. (2020b) Zhang, S.-N., Li, T., Lu, F., et al. 2020b, Science China Physics, Mechanics, and Astronomy, 63, 249502, doi: 10.1007/s11433-019-1432-6
  * Zhang et al. (2019) Zhang, Y., Ge, M., Song, L., et al. 2019, ApJ, 879, 61, doi: 10.3847/1538-4357/ab22b1
  * Zhang et al. (2022) Zhang, Y., Méndez, M., García, F., et al. 2022, MNRAS, 512, 2686, doi: 10.1093/mnras/stac690
  * Zhou et al. (2022) Zhou, D.-K., Zhang, S.-N., Song, L.-M., et al. 2022, MNRAS, 515, 1914, doi: 10.1093/mnras/stac1789