# A Phase-Resolved View to the Low-frequency Quasi-periodic Oscillations from
the Black Hole Binary MAXI J1820+070

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  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  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  A. Santangelo  Institut für
Astronomie und Astrophysik, Kepler Center for Astro and Particle Physics,
Eberhard Karls, Universität, Sand 1, D-72076 Tübingen, Germany  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  Y. Huang  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  Hong H. Liu
Center for Field Theory and Particle Physics and Department of Physics, Fudan
University, 200438 Shanghai, People’s Republic of China  P. Zhang  College of
Science, China Three Gorges University, Yichang 443002, China  W. Yu  Key
Laboratory of Particle Astrophysics, Institute of High Energy Physics, 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  J. Li  CAS Key Laboratory for Research in
Galaxies and Cosmology, Department of Astronomy, University of Science and
Technology of China, Hefei 230026, China  School of Astronomy and Space
Science, University of Science and Technology of China, Hefei 230026, China
Wen T. Ye  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  Pan P. Li  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  Zhuo L. Yu  Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 100049,
Beijing, China  Z. Yan  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  Yunnan Observatories, Chinese Academy of Sciences,
Kunming 650216, P.R. China

###### Abstract

Although low frequency quasi-periodic oscillations (LFQPOs) are commonly
detected in the X-ray light curves of accreting black hole X-ray binaries
(BHXRBs), their origin still remains elusive. In this study, we conduct phase-
resolved spectroscopy in a broad energy band for LFQPOs in MAXI J1820+070
during its 2018 outburst, utilizing Insight-HXMT observations. By employing
the Hilbert-Huang Transform method, we extract the intrinsic QPO variability,
and obtain the corresponding instantaneous amplitude, phase, and frequency
functions for each data point. With well-defined phases, we construct QPO
waveforms and phase-resolved spectra. By comparing the phase-folded waveform
with that obtained from the Fourier method, we find that phase-folding on the
phase of the QPO fundamental frequency leads to a slight reduction in the
contribution of the harmonic component. This suggests that the phase
difference between QPO harmonics exhibits time variability. Phase-resolved
spectral analysis reveals strong concurrent modulations of the spectral index
and flux across the bright hard state. The modulation of the spectral index
could potentially be explained by both the corona and jet precession models,
with the latter requiring efficient acceleration within the jet. Furthermore,
significant modulations in the reflection fraction are detected exclusively
during the later stages of the bright hard state. These findings provide
support for the geometric origin of LFQPOs and offer valuable insights into
the evolution of the accretion geometry during the outburst in MAXI J1820+070.

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

##  1 Introduction

Black hole X-ray binaries (BHXRBs) are commonly observed as transient sources,
with their primary feature being undergoing outbursts occasionally. Throughout
an outburst, the X-ray emission properties of the system evolves significantly
(Remillard & McClintock, 2006; Done et al., 2007). Most complete outbursts
typically exhibit four distinct canonical states: the low/hard state (LHS),
hard intermediate state (HIMS), soft intermediate state (SIMS) and high soft
state (HSS). Each state is characterized by distinct X-ray spectral and
variability properties (Belloni et al., 2005; Homan & Belloni, 2005). In
addition to X-ray emissions, BHXRBs are also characterized by radio/infrared
emissions, which are generally believed to be associated with relativistic
jets (Fender, 2001; Fender et al., 2004; Méndez et al., 2022, and references
therein). In the hard state, the radio emission with a flat spectrum is
interpreted as the self-absorbed synchrotron radiation from an optically
thick, steady and compact jet (Blandford & Königl, 1979; Fender, 2001). During
the hard-to-soft transition, the compact jet is gradually quenched, and the
radio emission is believed to originate from optically thin synchrotron
radiation produced by transient plasma clouds with relativistic velocities
(Fender et al., 2004).

The X-ray spectrum of a BHXRB is typically composed of multiple components.
First, there is a thermal component originating from a geometrically thin and
optically thick accretion disc (Shakura & Sunyaev, 1973; Lynden-Bell &
Pringle, 1974). Additionally, a non-thermal component arises from the Compton
up-scattering of disc photons in a hot corona region ($\sim 100$ keV) (Sunyaev
& Titarchuk, 1980; Titarchuk, 1994; Zdziarski et al., 1996, 1999; Życki et
al., 1999), possibly extending to the jet base (Markoff et al., 2005; You et
al., 2021). Furthermore, there is a reflection component resulting from a
fraction of the Comptonized photons illuminating the disc and being scattered
into the line of sight. The reflection spectrum exhibits prominent features
such as broad Fe K$\alpha$ lines at $\sim 6.4$ keV and a Compton hump in the
$\sim$20–30 keV energy range, etc. (Lightman & Rybicki, 1980; Fabian et al.,
1989; Dauser et al., 2010; García et al., 2014, 2015, and references therein).

In the time domain, low frequency quasi-periodic oscillations (LFQPOs, roughly
0.1–30 Hz, van der Klis, 1989) are usually seen in the light curves of BHXRBs,
characterized by a narrow peak with finite width in the power density spectra
(PDS). Based on the centroid frequency, quality factor and root-mean-square
(rms) amplitude, LFQPOs are often classified into three types: A, B and C
(Wijnands et al., 1999; Casella et al., 2005; Remillard & McClintock, 2006).
type-C QPOs are frequently observed in the LHS and HIMS, characterized by
strong amplitudes (fractional rms $\sim 10\%$) and the presence of flat-top
noise components in the PDS.

Over the past few decades, several models have been proposed to explain the
origin of LFQPOs in BHXRBs, based either on the geometric or the intrinsic
properties 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.. _Geometric_ models, on the
other hand, are mainly associated with the relativistic Lense-Thirring (L-T)
precession, which was initially proposed as the dynamic mechanism for QPOs by
Stella & Vietri (1998). As an extension of the relativistic precession model
(RPM; Stella et al., 1999), the L-T precession model proposed by Ingram et al.
(2009) assumes the entire hot flow precesses within the inner radius of the
truncated disc (Esin et al., 1997). Another _geometric_ model is the jet
precession model (Stevens & Uttley, 2016; de Ruiter et al., 2019; Ma et al.,
2021, 2023), supported by the recent simulation (Liska et al., 2018) and
observations of synchronous QPO signals in X-ray and optical/IR bands
(Kalamkar et al., 2016; Thomas et al., 2022).

Observational studies have revealed that the variability 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 prominent disc-like
component exists in the rms spectra (Sobolewska & Życki, 2006; Axelsson et
al., 2014; Axelsson & Done, 2016). Additionally, Ma et al. (2021) reported the
discovery of tpye-C QPOs above 200 keV, so the type-C QPO should be strongly
related to the Comptonized emission. Moreover, the inclination dependence of
amplitudes and time lags (see Motta et al., 2015; van den Eijnden et al.,
2017) and reflection variability extracted from the phase-resolved
spectroscopy (see Ingram & van der Klis, 2015; Ingram et al., 2016, 2017;
Nathan et al., 2022) provide further support for a geometrical origin of
LFQPOs, such as the L-T precession model. We refer readers to Ingram & Motta
(2019) for recent reviews of observations and theories of LFQPOs.

Phase-resolved analysis provides valuable insights into the properties of QPOs
by examining the phase-dependent behavior of the spectrum. In particular, the
L-T precession model predicts that reflection features, such as emission
lines, are modulated over the precession period (Ingram & Done, 2012; Ingram
et al., 2016; You et al., 2020). Therefore, conducting phase-resolved analysis
of QPOs offers a powerful diagnostic tool to distinguish between different
models. One possible method for generating phase-resolved spectra for QPOs is
the Hilbert-Huang Transform (HHT Huang et al., 1998). The HHT is a powerful
tool for analyzing phenomena with non-stationary periodicity and has been
successfully applied in astronomical researches. For instance, it has been
used to study the superorbital modulation of SMC X-1 (Hu et al., 2011), the
11-year sunspot variability (Barnhart & Eichinger, 2011), gravitational wave
signals (Camp et al., 2007; Hu et al., 2022), as well as QPOs in the active
galactic nucleus RE J1034+396 (Hu et al., 2014), and BHXRBs XTE J1550-564 (Su
et al., 2015) and MAXI J1820+070 (Yu et al., 2023). The HHT enables us to not
only trace the variation in frequency in the QPO but also to process phase-
resolved analyses even though the periodicity is unstable.

MAXI J1820+070 (ASASSN-18ey) is a galactic low-mass X-ray binary (LMXB)
discovered in optical by All-Sky Automated Survey for Super Novae (ASAS-SN,
Shappee et al., 2014) on 2018 March 3. Then it was discovered as an X-ray
transient by MAXI/GCS on 2018 March 11 (Denisenko, 2018; Kawamuro et al.,
2018; Tucker et al., 2018), located at R.A.=$18^{\rm h}20^{\rm m}21^{\rm s}.9$
and decl.=$+07^{\circ}11^{\prime}07^{\prime\prime}.3$ (Kawamuro et al., 2018).
So far, it is one of the brightest X-ray transients, having a soft X-ray flux
of $\sim 4$ Crab in 2–4 keV (Shidatsu et al., 2019) with a low column density
of $1.5\times 10^{21}\rm cm^{-2}$ (Uttley et al., 2018). Based on the radio
parallax measurement, Atri et al. (2020) determined a distance to the source
of $2.96\pm 0.33$ kpc. The VLBI observation provides an inclination angle of
$63\pm 3^{\circ}$ and a mass of $9.2\pm 1.3M_{\odot}$ (Atri et al., 2020).
With more accurate measurement of the mass ratio, Torres et al. (2020) gave a
mass of $8.48^{+0.79}_{-0.72}M_{\odot}$. Guan et al. (2021) and Zhao et al.
(2021) estimated the BH spin as $0.2^{+0.2}_{-0.3}$ and $0.14\pm 0.09$,
respectively, via modelling the continuum spectra. In this study, we
investigate LFQPOs in MAXI J1820+070 in details by utilizing the HHT method to
perform the phase-resolved spectral analysis. In Section 2, we describe the
Insight-HXMT observations and data reductions. In Section 3, we firstly
provide a concise introduction to the HHT method and its application in phase-
resolved spectral analysis of QPOs, and then present the corresponding
results. Subsequently, we discuss these results in Section 4 and summarize in
Section 5.

##  2 Observations and Data Reduction

In this study, we analyze data collected by the Hard X-ray Modulation
Telescope (Insight-HXMT) (Zhang et al., 2014, 2020b), which was launched on
June 15, 2017, as the first Chinese X-ray astronomy satellite. 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). These telescopes can operate in scanning and
pointing observational modes as well as in GRB mode. 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).

Insight-HXMT detected numerous type-C QPOs in MAXI J1820+070 during the hard
state of its 2018 outburst (see Ma et al., 2021; Yang et al., 2022). In this
work, we select 30 observations that span over the bright hard state, in which
the QPO frequency evolves from $\sim 0.04$ to $\sim 0.5$ Hz. In Figure 1, we
plot the Hardness-intensity diagram (HID) for MAXI J1820+070 during its 2018
outburst. Based on the observational dates and evolution of the count rate, we
categorize the selected 30 observations into four distinct epochs (plotted as
colored stars in Figure 1). Specifically, Epochs 1, 2, 3 and 4 correspond to
the time intervals MJD 58,204–MJD 58,209, MJD 58,220–MJD 58,225, MJD
58,228–MJD 58,239, and MJD 58,242–MJD 58,253, respectively. From Epochs 1 to
4, the count rate of LE evolves from $\sim 1000$ to $\sim 700$ cts s-1. As a
result, we combine more observations in Epochs 3 and 4 compared to Epochs 1
and 2, in order to achieve comparable count statistics within these four
epochs for subsequent QPO phase-resolved spectral analysis. Detailed
information regarding the aforementioned observations can be found in Table 1.
We extract the event data using Insight-HXMT Data Analysis Software v2.05,
along with the current calibration model
v2.06111http://hxmtweb.ihep.ac.cn/software.jhtml and the standard Insight-HXMT
Data Reduction Guide v2.05222http://hxmtweb.ihep.ac.cn/SoftDoc.jhtml, and
apply the following series of criteria recommended by the Insight-HXMT team:
(1) the elevation angle (ELV) is greater than $10^{\circ}$; (2) the geometric
cutoff rigidity (COR) is greater than 8 GeV; (3) the pointing position offset
is less than $0.04^{\circ}$; and (4) the good time intervals (GTIs) are at
least 300 s away from the South Atlantic Anomaly (SAA). The backgrounds are
produced 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 photon are corrected to the
barycenter of the solar system for their arrival times with the HXMTDAS tool
hxbary, and then evenly light curve with a bin size of 0.1 s for the
subsequent timing and spectral analyses of QPOs. The spectral analysis is
performed using the XSPEC v12.12.0 software package (Arnaud, 1996).

Figure 1: The hardness-intensity diagram (HID) of Insight-HXMT, where the hardness ratio is defined as the count ratio between the hard band (5–10 keV) and the soft band (1–5 keV), and the count rate is measured in the 1–10 keV energy range. The colored stars indicate the averaged results for each epoch. Table 1: Log of Insight-HXMT Observations Used in This Work # | ObsID | Start Time (MJD) | Exposure (s) | QPO Frequency (Hz) | Note
---|---|---|---|---|---
1 | P0114661006 | 58204.56 | 3325.7 | $0.046\pm 0.001$ | Epoch 1
2 | P0114661007 | 58205.55 | 2550.6 | $0.048\pm 0.003$ | Epoch 1
3 | P0114661008 | 58207.04 | 4649.3 | $0.051\pm 0.002$ | Epoch 1
4 | P0114661009 | 58207.83 | 3728.7 | $0.050\pm 0.003$ | Epoch 1
5 | P0114661010 | 58209.03 | 5356.6 | $0.055\pm 0.002$ | Epoch 1
6 | P0114661019 | 58220.82 | 9654.8 | $0.090\pm 0.004$ | Epoch 2
7 | P0114661020 | 58222.61 | 7601.0 | $0.103\pm 0.002$ | Epoch 2
8 | P0114661021 | 58223.20 | 2992.5 | $0.101\pm 0.003$ | Epoch 2
9 | P0114661022 | 58224.48 | 2573.5 | $0.110\pm 0.003$ | Epoch 2
10 | P0114661023 | 58225.48 | 4189.5 | $0.118\pm 0.004$ | Epoch 2
11 | P0114661026 | 58228.76 | 6492.8 | $0.135\pm 0.004$ | Epoch 3
12 | P0114661030 | 58229.89 | 2137.6 | $0.142\pm 0.005$ | Epoch 3
13 | P0114661031 | 58230.95 | 3247.9 | $0.155\pm 0.004$ | Epoch 3
14 | P0114661027 | 58231.94 | 2813.0 | $0.168\pm 0.005$ | Epoch 3
15 | P0114661028 | 58233.73 | 4639.4 | $0.173\pm 0.005$ | Epoch 3
16 | P0114661029 | 58235.53 | 4215.4 | $0.199\pm 0.006$ | Epoch 3
17 | P0114661032 | 58236.74 | 4029.9 | $0.259\pm 0.011$ | Epoch 3
18 | P0114661033 | 58238.04 | 1436.5 | $0.254\pm 0.007$ | Epoch 3
19 | P0114661034 | 58239.27 | 1460.4 | $0.278\pm 0.007$ | Epoch 3
20 | P0114661037 | 58242.22 | 1691.8 | $0.316\pm 0.009$ | Epoch 4
21 | P0114661038 | 58243.60 | 2871.8 | $0.329\pm 0.008$ | Epoch 4
22 | P0114661039 | 58244.61 | 1878.3 | $0.337\pm 0.008$ | Epoch 4
23 | P0114661040 | 58245.61 | 2719.2 | $0.409\pm 0.013$ | Epoch 4
24 | P0114661041 | 58246.60 | 5513.2 | $0.397\pm 0.007$ | Epoch 4
25 | P0114661042 | 58247.73 | 2632.4 | $0.461\pm 0.007$ | Epoch 4
26 | P0114661043 | 58248.72 | 2633.4 | $0.433\pm 0.011$ | Epoch 4
27 | P0114661044 | 58250.24 | 3950.1 | $0.426\pm 0.011$ | Epoch 4
28 | P0114661045 | 58251.03 | 2729.2 | $0.466\pm 0.010$ | Epoch 4
29 | P0114661046 | 58252.56 | 2394.0 | $0.462\pm 0.009$ | Epoch 4
30 | P0114661047 | 58253.15 | 1137.1 | $0.472\pm 0.011$ | Epoch 4

##  3 Analysis and Results

###  3.1 Extract Intrinsic QPO Light Curve Using Variational Mode
Decomposition

Phase-resolved spectral analysis has been extensively performed in studies of
eclipses and neutron star pulsations by folding the light curve. However, when
dealing with QPO signals in black holes, simply folding the light curve is not
appropriate due to the non-linear and non-deterministic evolution of their
phase over time (e.g. Morgan et al., 1997; Hu et al., 2014; Ingram & van der
Klis, 2015). The Hilbert-Huang transform, an adaptive data analysis method
developed by Huang et al. (1998), is a powerful tool to study signals with
non-stationary periodicity (Huang et al., 1998; Huang & Wu, 2008). This
approach consists of two main components: mode decomposition and Hilbert
spectral analysis (HSA). The mode decomposition aims to decomposing a time
series into several intrinsic mode functions (IMFs), while the HSA allows
obtaining both the frequency and phase functions of the desired IMFs, such as
the QPO component (e.g. Hu et al., 2014; Hsieh & Chou, 2020). The original
method for mode decomposition is the empirical mode decomposition (EMD)
proposed by (Huang et al., 1998), which directly operates in the temporal
space and adaptively decomposes the input signal into several IMFs. For a
detailed algorithmic description of EMD, we recommend referring to the review
by Huang & Wu (2008). The EMD is known to have limitations, such as the lack
of mathematical theory and sensitivity to noise and sampling. In contrast, the
variational mode decomposition (VMD) is a novel signal processing method that
surpasses traditional decomposition approaches of e.g. EMD (Dragomiretskiy &
Zosso, 2014). Compared to the traditional EMD technique, VMD theoretically
resolves the issue of mode mixing by decomposing the signal into a sum of IMFs
with analytically calculated limited center frequency and bandwidth. In the
VMD method, the IMFs are defined as amplitude-modulated-frequency-modulated
(AM-FM) signals, expressed as:

$u_{k}(t)=A_{k}(t)\cos{(\phi_{k}(t))},$ (1)

here the phase of the $k$-th mode, denoted as $\phi_{k}(t)$, is a non-
decreasing function, indicating that the instantaneous frequency,
$\omega_{k}(t)=\phi^{\prime}_{k}(t)\geq 0$. Additionally, the envelope
$A_{k}(t)$ is non-negative ($A_{k}(t)\geq 0$), while both the instantaneous
frequency $\omega_{k}(t)$ and envelope $A_{k}(t)$ exhibit much slower
variations compared to the phase $\phi_{k}(t)$ (Daubechies et al., 2011;
Gilles, 2013). Dragomiretskiy & Zosso (2014) demonstrated that a defined IMF
is a signal with limited bandwidth, which serves as the main assumption for
mode decomposition in the VMD. In other words, VMD assumes that each IMF is
primarily concentrated around a center frequency $\omega_{k}$. For a mode,
$u_{k}(t)$, with a mean frequency, $w_{k}$, its practical bandwidth increases
both with the maximum deviation $\Delta f$ of the instantaneous frequency from
its center, and with the rate of this excursion, $f_{\rm FM}$, according to
Carson’s rule (Carson, 1922). In addition to this, the bandwidth of the
envelope, $A_{k}(t)$, modulating the amplitude of the FM signal, itself given
by its highest frequency, $f_{\rm AM}$, broadens the spectrum even further.
Therefore, the total practical bandwidth of an IMF can be estimated as:
$BW=2(\Delta f+f_{\rm FM}+f_{\rm AM})$. To determine the bandwidth of each
mode, VMD utilizes several signal processing tools. These include applying the
Wiener filter for signal denoising, employing the Hilbert transform to
construct a single-sideband analytic signal, and performing frequency shifting
to baseband using complex harmonic mixing. We recommend to refer to
Dragomiretskiy & Zosso (2014) for detailed information regarding these signal
processing tools. The candidate mode is initially made single-sided, then
shifted to baseband, and its bandwidth is estimated using $H^{1}$ Gaussian
smoothness of the shifted signal, i.e. the $L^{2}$ norm of the gradient333In
Equation (2), $||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}$.:

$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},$ (2)

where $BW_{k}$ indicates the bandwidth of the mode, $\delta$ is the Dirac
function, and ‘$*$’ denotes convolution. The term enclosed in the square
brackets represents the analytic signal of $u_{k}(t)$, which will be
introduced in details in the subsequent text. For an input time series,
$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 objective
is to minimize the sum of bandwidths of the modes by solving the following
optimization problem:

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

where $K$ is the total number of modes. The first term in the objective
function promotes the compactness of the modes, while the second term measures
the reconstruction error that needs to be minimized. The parameter $\alpha$
serves as a weighting factor that balances the compactness and reconstruction
error. The VMD algorithm iteratively updates $\\{u_{k}\\}$ and
$\\{\omega_{k}\\}$ until convergence. In this process, $\omega_{k}$ is
directly calculated in the Fourier domain as:

$\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},$
(4)

where $\hat{u}_{k}(\omega)$ represents the Fourier transform of the mode
$u_{k}(t)$. It has been demonstrated that the solution of the optimization
problem can be considered as a generalization of the Wiener filtering into
adaptive, multiple bands. For more detailed description of the method, the
complete constrained variational optimization problem is available in
Dragomiretskiy & Zosso (2014). The choice of $K$ and $\alpha$ is crucial in
the VMD algorithm and determines its performance. Typically, higher values of
$K$ are associated with higher values of $\alpha$, indicating that as the
number of modes increases, the bandwidth of each mode decreases. In this
study, our ultimate goal is to perform the phase-resolved spectral analysis
for QPOs. Therefore, in the step of mode decomposition, our primary objective
is to extract the intrinsic light curve of the QPO fundamental as one IMF. To
accomplish this, we initially set $\alpha$ to a common value of 2000 and chose
2 modes ($K=2$). Subsequently, we calculate the power spectrum for each mode
and compare it with the power spectrum of the original light curve. If none of
the modes exhibit a PDS that matches the characteristic profile of the QPO, we
increase the value of $K$. This iterative process continues until we
successfully capture QPO’s fundamental and harmonic modes in an individual
mode. Finally, we adjust the value of $\alpha$ to match the profile of the PDS
of the QPO IMF with that of the QPO peak in the PDS of the original
lightcurve. In this work, 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. Figure 2 shows an example of a 100-s-long light
curve from LE (obsID P0114661038) with its corresponding IMFs. The light curve
exhibits a LFQPO with a frequency of $\sim$0.3 Hz. In this example, we set
$K=5$ and $\alpha=4000$ for the VMD algorithm. A $\sim 0.3$ Hz oscillation and
its corresponding harmonic component are identified as the second (plotted in
red) and third IMFs (plotted in blue), respectively. To provide a comparison
with the Fourier analysis, we present power density spectra of the original LE
light curve, as well as the second and third IMFs in the frequency range of
0.004 to 5 Hz in Figure 3.

Figure 2: Representative example of a 100-s-long lightcurve of LE in energy
range of 1–10 keV and its corresponding IMFs. In this example, we set $K=5$
and $\alpha=4000$ for the VMD algorithm, hence obtain five IMFs. The IMF C2
corresponds to the QPO signal at $\sim 0.3$ Hz, and the IMF C3 corresponds to
the second harmonic of the QPO. In this plot, the characteristic frequency
($f$) and the approximate bandwidth ($BW$) for each mode are also presented.
Figure 3: Power density spectra (PDS) of the original light curve (black
dots), IMF C2 (red line), and IMF C3 (blue line) plotted in Figure 2.

###  3.2 Hilbert Transform and QPO Waveform

For each IMF, the Hilbert transform can provide physically meaningful phase,
amplitude and frequency functions. The Hilbert transform of a time series
$f(t)$ is defined as

$\mathcal{H}(f(t))=\frac{1}{\pi}{\rm p.v.}\int{\frac{f(\tau)}{t-\tau}d\tau},$
(5)

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}(f(t))=A(t)e^{j\phi(t)},$ (6)

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}(f(t))\right]^{1/2}$ (7)

and

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

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

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

Considering the focus of the present study on investigating QPOs, we solely
apply the Hilbert spectral analysis to the QPO IMF (e.g. the second IMF
plotted in Figure 2). The left panel of Figure 4 illustrates the Hilbert
spectrum, represented by a color map showcasing the instantaneous frequency
and amplitude of the QPO component. The color depth represents the amplitude.
For comparison, we also employ the Weight Wavelet Z-transform (WWZ, Foster,
1996) technique for the time-frequency analysis. The right panel of Figure 4
displays the color contour of the WWZ power spectrum. In the WWZ analysis, we
adopt specific parameters, including a restricted frequency range of 0.08–1
Hz, a frequency step of 0.01 Hz, and a decay constant of $c=0.005$ (details on
these parameter settings can be found in Foster, 1996; Li et al., 2023). By
comparing the Hilbert spectrum with the WWZ power map, it is evident that the
QPO frequency fluctuates between 0.2 and 0.4 Hz over time. Importantly, the
main structures obtained from these two distinct methods are in agreement.
Both two methods show that, the $\sim$0.3 Hz oscillation is relatively weak in
time intervals of $\sim$ 18–25 s and $\sim$ 65–80 s. However, it becomes
stronger during time intervals of $\sim$ 40–60 s and $\sim$ 80–90 s, with
frequencies of $\sim 0.27$ and $\sim 0.35$ Hz, respectively (see Figure 4).

Figure 4: Left: Hilbert spectrum of $\sim 0.3$ Hz QPO IMF plotted in third
panel of Figure 2. The color on the z-axis represents the amplitude. Right:
two-dimensional contour map of the WWZ power spectrum. Both of the two
distinct methods present that the QPO frequency fluctuates between 0.2 and 0.4
Hz over time.

Once the instantaneous phases, $\phi(t)$, of the QPOs have been computed using
the Hilbert transform as described in Equation (8), the QPO waveforms can be
constructed by phase-folding the LE light curves (also see Hu et al., 2014;
Hsieh & Chou, 2020). To generate the QPO waveform, we divide the phases into
30 bins. Figure 5 illustrates the results for the four epochs. The QPO
waveform exhibits a distinctive non-sinusoidal nature, characterized by a
gradual rise of $\sim$0.7 cycles, followed by a rapid drop of $\sim$0.3
cycles. This non-sinusoidal behavior is associated with the presence of the
second harmonic in the PDS (see Figure 3).

Figure 5: Constructed QPO waveforms of four epochs by phase-folding the
lightcurve from LE in the energy range of 1–10 keV.

###  3.3 Phase-resolved Spectra

With the well-defined phase, we construct Good Time Intervals (GTIs)
corresponding to 10 distinct phase intervals. These GTIs enable us to extract
spectra for subsequent phase-resolved spectral analysis in the four epochs.
For the spectral analysis, energy bands of 2–10 keV for LE, 10–30 keV for ME,
and 30–120 keV for HE, respectively, are adopted. Due to calibration
uncertainties at 22 keV in the ME, we exclude the 20–23 keV range during
spectral analysis for this instrument. To account for the current accuracy of
the instrument calibration, we introduce systematic errors of 0.5%, 0.5%, and
1% to the energy spectra for LE, ME, and HE, respectively. We employ the same
spectral model as in You et al. (2021), which is $constant\times
Tbabs\times(diskbb+relxillCp+xillverCp)$, to fit these spectra at different
phases. The constant factor is utilized to accommodate flux calibration
discrepancies among the three distinct instruments. In phase-averaged spectral
fittings, we set the constant factor of LE to 1 and treat those of ME and HE
as free parameters, allowing them to vary during the fittings. The $Tbabs$
accounts for interstellar absorption, with the equivalent hydrogen column,
$N_{\rm H}$, fixed at $0.15\times 10^{22}$ cm-2 (Uttley et al., 2018). In the
non-relativistic multi-color blackbody model $diskbb$, two crucial parameters
are estimated, namely the temperature ($T_{\rm in}$) and the inner radius
($R_{\rm in}$) of the accretion disk. The latter one is linked to the
normalization of the $diskbb$ model follows as $N_{\rm disk}=(R_{\rm
in}/D_{10})^{2}\cos i$, where $N_{\rm disk}$ is the normalization, $R_{\rm
in}$ is in units of km, $D_{10}$ is the distance to the source in units of 10
kpc and $i$ is the inclination. By taking the black hole mass as $8.48\
M_{\odot}$ (Torres et al., 2020), the distance and inclination as 2.96 kpc and
63∘ , respectively (Atri et al., 2020), we can obtain $N_{\rm disk}\approx
816(R_{\rm in}/R_{\rm g})^{2}$, where $R_{\rm g}=GM/c^{2}$ is the
gravitational radius. For the relativistic reflection model $relxillCp$, we
assume the canonical power-law emissivity profile, $\epsilon=r^{-3}$, for the
disk (Fabian et al., 1989). Referring to previous studies, the disk
inclination ($i$), spin parameter of the black hole ($a*$) are set to
$63^{\circ}$, 0.14, respectively (Atri et al., 2020; Zhao et al., 2021). The
parameter $R_{\rm out}$ is the outer radius of the accretion disk, which turns
out to not be sensitive to the overall fitting and hence is frozen at the
maximum value, 1000$R_{\rm g}$. Since $R_{\rm in}$ can also be a free
parameter in $relxillCp$, it is self-consistent to link it to $N_{\rm disk}$
in $diskbb$ as $N_{\rm disk}=816R_{\rm in}^{2}$ (see also Gao et al., 2023),
where $R_{\rm in}$ is in units of $R_{\rm g}$. However, it is important to
emphasize that, for the sake of simplicity, we have disregarded any
relativistic effects when establishing the connection between the two
parameters. The electron density ($N_{\rm e}$) is set to the default value of
1015 cm-3 since our energy coverage (2–120 keV) does not allow for the
constraint of soft excess emission below 2 keV (García et al., 2016). In the
spectral fitting, we also apply the non-relativistic reflection model
$xillverCp$ to account for the distant reflector with the narrow/low-
ionization iron line (e.g., García et al., 2015). The parameters in
$xillverCp$ are tied to the corresponding ones in $relxillCp$ except for the
ionization parameter and the normalization. Since the spectral fits are
insensitive to the ionization parameter, we fix the ionization parameter at a
low value with $\log\xi=1.0$ (see also You et al., 2021). Hence, the
normalization of $xillverCp$ is the only free parameter in the model. In our
phase-resolved spectral fittings, we fix constant values, the iron abundance
$A_{\rm Fe}$ and the ionization parameter $\log\xi$ to the same values
obtained from the phase-averaged spectral fittings. After fitting the energy
spectra with the model, we obtain a satisfactory reduced $\chi^{2}$ ($\sim
1$). To compute the uncertainties at the 90% confidence level, we employ the
Markov Chain Monte Carlo (MCMC) method using the Goodman-Weare algorithm with
8 walkers and a total length of 40,000 (Goodman & Weare, 2010), and the
initial 2000 elements are discarded as “burn-in” period. We find that the
autocorrelation length is typically one thousand elements, so the net number
of independent samples in the parameter space we have is order of $10^{4}$. In
order to further test the convergence of the MCMC chain, we compare the one-
and two-dimensional projections of the posterior distributions for each
parameters from the first and second halves of the chain (see also Liu et al.,
2021), and we find no significant differences. In Appendix A, we present the
contour maps and probability distributions for each free parameter.

Figure 6: Phase dependence of the seven free spectral parameters, including
the temperature of the inner radius of the disk ($T_{\rm in}$), inner radius
($R_{\rm in}$), spectral index ($\Gamma$), electron temperature ($kT_{\rm
e}$), reflection fraction ($R_{\rm f}$), and the normalization of $relxillCp$
(Norm1) and $xillverCp$ (Norm2), across the four epochs. Additionally, we
include the QPO waveform in each panel as the gray dashed line.

Figure 6 displays the phase dependence of the seven free parameters, including
the temperature of the inner radius of the disk ($T_{\rm in}$), inner radius
($R_{\rm in}$), spectral index ($\Gamma$), electron temperature ($kT_{\rm
e}$), reflection fraction ($R_{\rm f}$), and the normalization of $relxillCp$
(Norm$1$) and $xillverCp$ (Norm$2$), across the four epochs. Additionally, we
include the QPO waveform in each panel as the gray dashed line. From the
phase-resolved spectral analysis, we observe no significant modulations in
$T_{\rm in}$ and $R_{\rm in}$ as the flux varies on the QPO period. However,
notable modulations are observed in the spectral index and the normalization
of $relxillCp$ throughout all four epochs. Moreover, these two parameters
exhibit the same modulation phase as the flux. Additionally, significant phase
modulations are observed for parameters $R_{\rm f}$ and $kT_{\rm e}$ in Epochs
3 and 4, but not 1 and 2. For the modulation of $R_{\rm f}$ in Epochs 3 and 4,
it seems to lead the flux by $\sim$0.2 cycles. On the other hand, the
normalization of $xillverCp$ evolves from strong phase modulation in Epochs 1
and 2 to relatively marginal modulation in Epochs 3 and 4.

##  4 Discussion

In this work, we have analyzed observations of MAXI J1820+070 from Insight-
HXMT and extracted the intrinsic QPO variability by utilizing the VMD
algorithm to the observed light curves. Through the application of the Hilbert
transform to the intrinsic QPO variability, the instantaneous amplitude, phase
and frequency functions are obtained for each data point. The HSA reveals that
the QPO exhibits variations in its frequency and amplitude. Via phase-revolved
spectral analyses we find strong and persistent phase modulation of the
spectral index, and moderate evolution of phase modulation in the reflection
fraction. These findings provide valuable insights into the origins of LFQPOs
and the evolution of accretion geometry during outbursts in BHXRB systems.

###  4.1 QPO Waveform

In Section 3, the QPO waveforms are constructed through the phase-folding of
light curves (see Figure 5). Despite that the non-sinusoidal nature,
characterized by a distinct “slow rise and fast decay” feature, is observed in
the QPO waveforms across the four epochs, the contribution of the harmonic
component in these waveforms is relatively weaker than that observed in PDS.
The phase difference between QPO fundamental and harmonic modes has been
demonstrated to be unstable in time (Ingram & van der Klis, 2015; de Ruiter et
al., 2019). This indicates that the phase-folding on the phase of the
fundamental mode could suppress the contribution of the harmonic component in
the derived waveform. To demonstrate this suppressing effect in the phase-
folding method, we compare the phase-folded waveform with that obtain from the
Fourier method. To do this, we conduct a similar analysis as Ingram & van der
Klis (2015) to reconstruct the QPO waveform for Epoch 4. In order to measure
the phase of the QPO fundamental and harmonic components, we first perform the
Fourier transform of many short segments of each light curve within Epoch 4.
The phases of the fundamental and harmonic components for a single segment are
calculated as the argument of the Fourier transform for that segment at the
frequency bins containing the QPO’s fundamental and harmonic components,
respectively. Specifically, for the harmonic component, we utilize the
frequency bin that is precisely twice the frequency of the bin used for the
fundamental component. The phase difference for each segment is calculated
using Equation (3) of Ingram & van der Klis (2015). The distribution of the
phase difference ($\psi/\pi$) is presented in the left panel of Figure 7,
where the results obtained from the Fourier method is plotted in red. On the
other hand, the HHT method enables us to obtain the instantaneous phase for
each IMF. Thus, we calculate the phase for the IMFs that corresponds to the
QPO’s fundamental and harmonic components, respectively, and then determine
the phase difference for each time bin ($dt=0.1$ s). The resulted $\psi/\pi$
distribution derived using the HHT method is plotted in black in the left
panel of Figure 7. As one can see, both methods yield similar $\psi/\pi$
distributions, but the HHT method produces a more concentrated distribution.
The mean values of $\psi/\pi$ from the Fourier and HHT methods are found to be
$0.191\pm 0.041$ and $0.229\pm 0.002$, respectively. We propose that the HHT
method could serve as an effective and novel approach to determine the phase
difference between QPO harmonics. Based on the RMS of the fundamental and
harmonic components, $\psi/\pi$, as well as the Equation (6) of Ingram & van
der Klis (2015), the QPO waveform can be reconstructed. The waveforms obtained
from both the phase-folding and Fourier methods are presented in the right
panel of Figure 7. Although both waveforms exhibit a non-sinusoidal nature,
characterized by a “slow rise and fast decay” feature, the contribution of the
harmonic component in the waveform derived from the HHT method is indeed
suppressed. It is worth noting that this suppressing effect may also introduce
slight differences in the modulations of spectral parameters between the HHT
and Fourier methods.

Figure 7: Left: Histogram of measured values of the phase difference between
harmonics obtained using the HHT (black) and Fourier (red) methods. To show a
clear peak, we repeat the data from 0 to $\pi$. Right: Reconstructed QPO
waveform derived the HHT phase-folding (black) and Fourier (red) methods.

de Ruiter et al. (2019) measured the phase difference for many observations of
QPOs in the RXTE archive. In most high-inclination systems, $\psi/\pi$ between
harmonics is distributed around $\sim 0$ to $\sim 0.4$. This suggests that the
majority of QPO waveforms observed in these systems exhibit the characteristic
of “slow rise and fast decay”, even though the evolution of $\psi$ may
slightly alter the waveform. Considering that the inclination of MAXI
J1820+070 is $\sim 63^{\circ}$, the QPO waveforms observed in these four
epochs, which display the “slow rise and fast decay” feature, are consistent
with the behavior observed in most high-inclination systems from RXTE
observations. In Epoch 4, we observe $\psi/\pi\sim 0.2$ for QPOs at $\sim 0.4$
Hz, which is slightly smaller than the values observed in other high-
inclination systems at similar frequencies. However, in order to make a
comprehensive comparison of the evolution of $\psi$ in MAXI J1820+070 with
that in other sources, it is necessary to have more QPO samples at higher QPO
frequencies ($\sim 1-5$ Hz) which were not detected by Insight-HXMT.

###  4.2 Modulation of the Spectrum

As shown in Figure 6, the spectral index ($\Gamma$) and the normalization of
$relxillCp$ (Norm1) exhibit significant modulations, which closely follow the
flux modulation. In contrast, the modulations of parameters in $diskbb$ are
marginal. This suggests that the dominant contribution to the QPOs comes from
the Comptonized component, supporting models such as the corona’s L-T
precession model (Ingram et al., 2009) and the time-dependent Comptonization
model (Karpouzas et al., 2020; Bellavita et al., 2022). Since a jet could also
act like a corona which produces Comptonized photons (Markoff et al., 2005),
the jet precession model is also a plausible explanation for the findings
(Stevens & Uttley, 2016; Liska et al., 2018; Ma et al., 2021, 2023).

####  4.2.1 Corona Precession Model

If the QPOs are produced by the corona’s Lense-Thirring precession, the X-ray
QPO variability is due to the geometric wobble of the corona, which changes
the projection area of the corona with respect to the observer to modulate the
observed X-ray flux (Ingram et al., 2009, 2015; You et al., 2018, 2020; Shui
et al., 2023).

Figure 8: Phase dependence of the optical depth across the four epochs. The
colored scatters represent the calculated optical depth ($\tau$) from Equation
10, using the best-fit parameters of $\Gamma$ and $kT_{\rm e}$ from the
spectral fitting. The black dashed lines represent the calculated optical
depth based on Equation 11 and 12.

Since the geometric shape of a corona might not be a sphere (e.g., a crushed
sphere, Ingram & Done, 2012; Shui et al., 2023), the optical depth of the
corona respect to the observer could vary over the precession period. The
illuminating continuum in $relxillCp$ is $nthComp$ (Życki et al., 1999),
wherein the spectral shape is determined by the combination of the electron
temperature ($kT_{\rm e}$) and the optical depth ($\tau$). Therefore, the
modulations of $\Gamma$ could be attributed to the modulations of $kT_{\rm e}$
and $\tau$. For a given $kT_{\rm e}$ and spectral index ($\Gamma$), we can
estimate the corona’s optical depth by

$\tau=\sqrt{\frac{9}{4}+\frac{m_{\rm e}c^{2}}{kT_{\rm
e}}\frac{3}{(\Gamma-1)(\Gamma+2)}}-\frac{3}{2},$ (10)

where $m_{\rm e}$ is the electron rest mass, $c$ is the speed of light
(Zdziarski et al., 1996). We utilize Equation (10) to calculate the optical
depth at different phases, using the best-fit parameters of $\Gamma$ and
$kT_{\rm e}$ from the spectral fittings. The modulations of the optical depth
for four epochs are plotted in Figure 8 as colored scatters. For all four
epochs, the optical depth exhibits contrasting modulations with respect to the
flux. This phenomenon can be qualitatively explained by the corona’s Lense-
Thirring precession model. When the corona is oriented face-on towards the
observer, the projected area and consequently the flux of the corona appear to
be larger compared to the edge-on orientation. However, the optical depth of
the corona modulates in the opposite manner due to the self-occultation effect
(You et al., 2018). Quantitatively, according to simplified calculations by
Shui et al. (2023), if the corona’s shape resembles a crushed sphere, there
might exist an anti-correlation between the flux and corona’s optical depth.
Equation (9) of Shui et al. (2023) describes the modulation of the corona’s
optical depth as the corona precesses:

$\tau=\tau_{0}\left\\{\left(h/r\right)^{2}+\left[1-\left(h/r\right)^{2}\right]\cos^{2}{\theta}\right\\}^{-1/2},$
(11)

where $\tau_{0}$ is the minimum optical depth of the corona, i.e. viewed from
the corona’s normal, and $h/r$ is the scale height which describes the shape
of the corona. When $h/r=1$, the corona’s shape is a sphere. $\theta$ is the
included angle between the line of sight and the corona’s normal, and it
varies over the precession period. The cosine of $\theta$ can be expressed
using the following equation:

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

where $i$ is the inclination angle between the line of sight of the observer
and the black hole spin axis, $\Phi$ is the azimuth angle of the observer,
$\gamma$ is the included angle between the corona’s normal and the black hole
spin axis and $\omega$ is the precession phase angle, which varies from 0 to
2$\pi$ during the precession period. We set $i=63^{\circ}$, $\Phi=180^{\circ}$
(resulting in the peak flux occurring at 0.5 QPO cycles) and
$\gamma=5^{\circ}$ to calculate modulations of the optical depth and plot the
results in Figure 8 as dashed black lines. It is evident that the theoretical
curves closely match the data points calculated using Equation (10) across all
four epochs. For the four distinct epochs, we adjust the values of $\tau_{0}$
and $h/r$ to make the calculated modulation curves comparable to the data
points. However, we note that the initial parameter values solely influence
the amplitude and mean values, without affecting the overall trend of the
modulations in $\tau$. It is important to note that the simplified model
mentioned above disregards the Doppler boosting effect. In reality, the
modulation of the X-ray flux is primarily influenced by a combination of the
projected area and Doppler boosting effects. However, maximal Doppler boosting
occurs when the flow is viewed most edge-on, since blue shifts from the
approaching material dominate over red shifts from receding material. On the
other hand, the projected area is largest when the flow is viewed most face-
on. Therefore, the X-ray flux as a function of precession phase represents the
competition between these two effects (Veledina et al., 2013; Ingram et al.,
2015). In Figure 8, we can observe a good alignment between the data points
and the theoretical curves calculated using Equation (11), which considers
only the projected area effect. This suggests that the modulation of the X-ray
flux across the four epochs is likely dominated by the projected area effect.

As shown in Figure 6, we detect significant modulations of reflection fraction
($R_{\rm f}$) in Epochs 3 and 4. Reflection fraction is defined as the ratio
of Comptonized emission intensity that illuminates the disk to the Comptonized
emission intensity that directly reaches the observer. This parameter is tied
to the intrinsic accretion geometry (Dauser et al., 2016). Therefore, the
observed modulation in $R_{\rm f}$ reflects that the accretion geometry keeps
changing over the cycle and thus the solid angle of the emitter as seen by the
reflector and/or the solid angle of the reflector as seen by the observer will
also change accordingly. Consequently, the observed modulations in reflection
fraction provide strong support for the geometric origin of the QPO. For
instance, within the framework of the corona’s L-T precession model proposed
by Ingram et al. (2009), it is extensively anticipated that the reflected
emission would display variations throughout the QPO period (Ingram & Done,
2012; Ingram & van der Klis, 2015; Ingram et al., 2016; You et al., 2018,
2020; Nathan et al., 2022), subsequently leading to corresponding modulations
in the derived reflection fraction (see You et al., 2020).

####  4.2.2 Jet Precession Model

Another geometric model worth considering is the jet precession model proposed
by Ma et al. (2021, 2023). In this model, the main picture is that, the
emission on-axis is significantly stronger than that off-axis due to the
Doppler boosting, since typical speeds at the base of the jet can already be
relativistic (McKinney, 2006). As a result, the observed flux from the
precessing jet can be modulated. In this case, the speed of the jet material
can effectively impact the QPO variability. Specifically, the modulation
amplitude increases with higher jet speeds. This is because larger speed leads
to greater differences in emission strength between the on-axis and off-axis
cases.

Ma et al. (2021), Zhou et al. (2022), Ma et al. (2023) and Yu et al. (2023)
have recently found that the soft energy photons intrinsically lag to the hard
ones for the QPO signals in MAXI J1820+070. In the context of the jet
precession model, these soft lags can be attributed to the differences among
emitting regions for the emitted photons in different energy bands within the
jet, as well as the geometric twist of the precessing jet itself. Typically,
the harder photons are emitted at the lower part of the jet, while softer
photons originate as the jet evolves further out. As the jet undergoes
precession, different regions of the jet experience different precession
phases, resulting in the observed soft lags. This original idea assumes that
the speeds of the jet materials in different regions are the same. However,
the jet at the base is at rest and from there on the particles are very
efficiently accelerated (McKinney, 2006). This means the jet does not have a
uniform velocity (see also Dauser et al., 2013). If the QPO is indeed produced
by the jet precession, the observed similarity in modulated phases between
$\Gamma$ and the total flux suggests that the jet speed increases from lower
to higher heights. As the speeds vary across different heights (corresponding
to different energy bands of the emission), the modulation amplitude at
specific phases differs for each energy band. Specifically, lower energy
photons are emitted from higher-speed ejecta, resulting in the normalized flux
of lower energies being larger than that of higher energies during the peak
phases. During the trough phases, the situation is reversed, where the
normalized flux of higher energies becomes larger than that of lower energies.
Therefore, the hardness should vary over the precession cycle. In Figure 9, we
plot the QPO waveforms for three instruments. To plot this figure, we perform
the HHT to the light curves from each instrument and phase-fold them, then
modify the phase differences among the three instruments. The plotted QPO
waveforms are normalized by the averaged fluxes. Additionally, we plot the
hardness, calculated as the ratio of the count rate in the HE band to that in
the LE band within each phase intervals, in Figure 9 as a black line. Clearly,
the hardness exhibits an inverse modulation with respect to the flux, and the
variability amplitude shows an increasing trend with energy, which is
consistent with RMS spectrum presented by Ma et al. (2021). However, it is
important to note that the maximum soft lags (photons at $\sim 100$ keV with
respect to that at $\sim 2$ keV) are $\sim$0.5 s at the QPO frequency of
$\sim$0.04 Hz (Ma et al., 2021), indicating a phase lag of $<0.05$ QPO cycles.
Consequently, the phase lags between different energy bands are not expected
to have a significant influence on the phase-resolved spectral analysis.
Therefore, the spectra obtained at the peak and trough phases in Section 3
allow us to establish the relationship between the speed of the ejecta and the
energy of the corresponding emitted photons. For a continuous jet, the
observed flux from the jet can be written by

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

where $F_{\rm obs}$ is the observed spectrum, $F_{\rm int}$ is the intrinsic
spectrum of the jet, $\Gamma$ is the spectral index 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}},$ (14)

where $\beta=v/c$ is the jet speed in units of the speed of light, and
$\theta$ is the included angle between the jet axis and the line of sight of
the observer, which is calculated by Equation (12). Therefore, the ratio of
the trough spectrum to peak spectrum can be expressed as

$\frac{F_{\rm t}}{F_{\rm p}}=\left[\frac{1-\beta\cos\theta_{\rm
p}}{1-\beta\cos\theta_{\rm t}}\right]^{\Gamma+2},$ (15)

where $F_{\rm p}$ and $F_{\rm t}$ are the spectra at the peak and trough
phases, respectively, and $\theta_{\rm p}$ and $\theta_{\rm t}$ are the
included angles at the peak and trough phases, respectively. In the case of
$\Phi=180^{\circ}$, $\theta_{\rm t}=\theta_{\omega=0}$ and $\theta_{\rm
p}=\theta_{\omega=\pi}$, respectively. We define the ratio between the trough
spectrum and peak spectrum as $\mathcal{R}=F_{\rm t}/F_{\rm p}$, which is a
function of energy, and plot the calculated results for four epochs in the top
panels of Figure 10. It can be observed that $\mathcal{R}$ exhibits an
increasing trend with energy across four epochs, indicating that the trough
spectra is harder than the peak spectra. From Equation (15), it can be deduced
that

$\beta(E)=\left[\mathcal{R}^{1/(\Gamma+2)}-1\right]\cdot\left[\mathcal{R}^{1/(\Gamma+2)}\cos\theta_{\rm
t}-\cos\theta_{\rm p}\right]^{-1},$ (16)

where $E$ is the energy of emitted photons. Here, we set $\Gamma$ as the value
obtained from the spectral fitting of the phase-averaged spectra. The bottom
panels of Figure 10 displays the energy dependence of $\beta$. The obtained
anti-correlation between $\beta$ and the energy indicates as the jet material
is ejected away from the black hole and the energy of the emitted photons
decreases, while the speed of the jet material increases. We also observe a
slight decrease in the computed jet speed in each energy band from Epochs 1 to
4. Additionally, Ma et al. (2021) presented a decreasing trend of QPO RMS
during the same stage. These findings collectively suggest that the outflowing
materials decelerate during the bright hard state according to the jet
precession model. However, the evolution of $R_{\rm f}$ indicates an
acceleration of the outflowing plasma (You et al., 2021). This apparent
contradiction may arise from a potential oversimplification within the jet
precession model, specifically in maintaining a constant included angle
between the jet and black hole spin axes ($\gamma$). In fact, the coupling
between the jet and inner disk (see Ma et al., 2023) could entail a decrease
in $\gamma$ as the inner radius of the disk moves inwards. This behavior
arises due to the tendency of the disk to align with the black hole spin,
driven by the Bardeen-Petterson (BP) effect (Bardeen & Petterson, 1975).
Consequently, the precession of the jet may occur with a reduced $\gamma$.
Given that the QPO RMS is influenced by both the jet speed and $\gamma$, the
combined impact of jet acceleration and the $\gamma$ reduction may account for
the observed declining trend in QPO RMS as well as the computed jet speed
using Eq. (16). It is important to note that the derived trend of $\beta$-$E$
correlation within each epoch remains unaffected by the assumption of a
constant $\gamma$.

Figure 9: Constructed QPO waveforms of Epoch 1 by phase-folding the lightcurve
from LE, ME and HE in the energy ranges of 1–10 keV (orange), 10–30 keV (blue)
and 30–120 keV (pink), respectively. The black line represents the hardness
defined as ratio between constructed QPO waveforms of HE and LE.

Figure 10: Top: the ratio between the spectra of QPO trough (0–0.1 cycles) and
peak (0.4–0.5 cycles) phases across the four epochs. Bottom: the energy
dependence of the jet speed ($\beta$) across the four epochs, where values of
$\beta$ are calculated by Equation 16. The gray scatters represent the
calculation from the original spectra, while colored scatters represent
rebinned values.

Using Insight-HXMT and NICER observations, Ma et al. (2023) presented the
energy dependence of the LFQPO characteristics and phase lags from 0.2 to 200
keV, which could be explained by the picture that the jet and inner disk ring
precessing together. In this model, the inner region of the jet and the
innermost part of the accretion disk are somehow coupled, resulting in a
scenario where the precession of the disk induces a corresponding precession
of the jet. Crucially, the geometrical alignment between the jet and the inner
disk remains constant throughout the precession cycle. In this case, if the
height of the X-ray emission region within the jet is sufficiently low,
thereby solely illuminating the precessing inner disk, the reflection fraction
($R_{\rm f}$) should in principle remain constant over the QPO cycle. However,
the inclination of the disk, denoted as $i$, is expected to undergo
modulations that contrast with the flux variations. Taking this into account,
we set $R_{\rm f}$ to the best-fit values acquired from the phase-averaged
spectral fittings, while simultaneously allowing the inclination angle, $i$,
to be a free parameter in the phase-resolved spectral fittings. Following the
adjustment to spectral fittings, we also obtain a reasonable value of
$\chi^{2}$ ($\sim$1), and compute the uncertainties of $i$ at the 90%
confidence using the MCMC technique. The behavior of $i$ across the QPO cycle
for the four epochs under investigation is depicted in Figure 11. Clearly, the
inclination exhibits notable modulations and shows large phase lags relative
to the flux in Epochs 3 and 4, whereas positive correlations between $i$ and
flux are observed in Epochs 1 and 2. This modulation pattern of $i$ during
Epochs 3 and 4 lends support to the aforementioned jet-disk coupling
precession model. On another note, our spectral fitting model employs
$xillverCp$ to accommodate the contribution arising from the distant reflector
component. Intriguingly, the modulation of $xillverCp$ normalization seems to
be more conspicuous during Epochs 1 and 2 in contrast to Epochs 3 and 4 (see
Figure 6). It is challenging to explain these distinctions across four epochs.
One plausible interpretation is that the underlying phase-averaged accretion
geometry undergoes changes during the bright hard state.

From a geometric perspective, a higher vertically extended jet would naturally
have a larger solid angle with respect to the accretion disk, thereby
efficiently illuminate the outer disk. As a result, both the inner precessing
and the outer non-precessing regions of the disk could effectively receive
incident flux originating from the vertically extended jet. However, it is
important to note that modulations in the reflection component from these two
regions are in principle different, owing to the fact that the inner disk
undergoes precession to align with the jet, while the outer disk remains
stable. Consequently, the resulting modulations in reflection fraction may
potentially undergo some degree of smearing out. Additionally, the reflection
model $relxillCp$ assumes a symmetric geometry (see e.g. Dauser et al., 2013),
while the illumination pattern in the above precession scenario is asymmetric.
Therefore, the positive correlation between $i$ and flux over the QPO cycle,
as shown in Epochs 1 and 2, might be attributed to unclear parameter couplings
within the current reflection models. The resolution of this matter warrants a
detailed investigation through further simulations.

As mentioned above, if the height of the X-ray emission region within the jet
is small enough to illuminate the inner precessing disk solely, modulations of
$R_{\rm f}$ can be readily detected. Because the inner precessing disk is the
primary reflector, leading to a more singular modulation pattern of reflection
component. Furthermore, the observed substantial phase differences between the
inclination angle ($i$) and the flux, align with the fundamental concept of
the jet-disk coupling precession model. We propose that the picture of the
decreasing jet emission height is consistent with the scenario of a
contracting Comptonization region during the bright hard state (Kara et al.,
2019; You et al., 2021; Wang et al., 2022). It is important to acknowledge
that the jet-disk (where the disk corresponds to the geometrically thin and
optically thick reflector) coupled precession model has not been supported by
recent general relativistic magneto-hydrodynamic (GRMHD) simulations. For
instance, Liska et al. (2018) demonstrated that the thick accretion flow and
jets precess together, owing to the combination of Lense–Thirring effects and
pressure or magnetic torques from the inflow/outflow system. However, their
simulation did not include the optically thick disk. More recently, a
simulation conducted by Bollimpalli et al. (2023) included both thin and thick
accretion disks, revealing that the thick disk undergoes precession while the
thin disk remains relatively stable. From our results, it remains uncertain
whether the thin disk can undergo precession. Therefore, further investigation
using numerical simulations and observations is necessary to explore this
aspect in the future.

Regarding the modulation of the distant reflector component, it is not at odd
that modulation strength decreases when the height of the jet becomes
significantly smaller in the jet-disk coupling precession model. However, it
is crucial to determine whether or not the modulation of Norm2 is physical.
Specifically, if the light-crossing time is longer than the QPO period, the
modulation of the distant reflection component will be washed out. For the
observed QPOs in Epoch 1, with a period of $\sim 20$ s, the corresponding
light-crossing distance is estimated to be $\sim 6\times 10^{6}$ km
(equivalent to $\sim 4.8\times 10^{5}R_{\rm g}$, assuming $M_{\rm
BH}=8.5M_{\odot}$). Based on the well determined orbital parameters of the
system, the outer disc radius is constrained to be between $\sim 6.3\times
10^{4}$ and $\sim 2.5\times 10^{5}R_{\rm g}$ (Koljonen et al., 2023).
Therefore, if the reflector is located within the distant accretion disk, it
is possible to observe modulation in the distant reflection component.
However, no significant phase lags between the total flux and this component,
as shown in Figure 6, may suggest that the modulation of Norm2 in Epochs 1 and
2 is not likely from reflection from the distant disk.

Figure 11: Phase dependence of the inclination across the four epochs. To
obtain the modulations in the inclination ($i$), we set reflection fraction to
the best-fit values acquired from the phase-averaged spectral fittings. The
uncertainties are reported at the 68% confidence using the MCMC technique.
Additionally, we include the QPO waveform in each panel as the gray dashed
line.

##  5 Conclusion

With the VMD algorithm and Hilbert transform applied to the LFQPOs of MAXI
J1820+070 observed by Insight-HXMT, the QPO waveforms are constructed. The
phase-resolved energy spectra show strong phase modulation of the spectral
index concurrent to the source flux. Similar phase modulations are found for
the reflection fraction and evolve with outburst. These findings provide new
insight to understanding the QPO in aspects of origination and evolution
within outburst.

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 MCMC Parameter Probability Distributions

In this section, we show an example (the QPO trough phase in Epoch 4) of the
MCMC analysis for the spectral fitting model presented in Section 3. We use
the Goodman-Weare algorithm with 8 walkers and a total length of 40,000 to
perform the MCMC analysis, and the initial 2000 elements are discarded as
“burn-in” period during which the chain reaches its stationary state. We find
that the autocorrelation length is typically one thousand elements, so the net
number of independent samples in the parameter space we have is order of
$10^{4}$. Furthermore, we compare the one- and two-dimensional projections of
the posterior distributions for each parameters from the first and second
halves of the chain to test the convergence, and we find no significant
differences (see Figure 12). The contour maps and probability distributions
are plotted using the corner package (Foreman-Mackey, 2016).

Figure 12: One- and two-dimensional projections of the posterior probability
distributions, and the 0.16, 0.5 and 0.84 quantile contours derived from the
MCMC analysis for each free spectral parameters. To test the convergence, we
compare the one- and two-dimensional projections of the posterior
distributions from the first (blue) and second (orange) halves of the chain,
and we find no large differences. This illustration corresponds to the
spectral fitting of the QPO trough phase in Epoch 4.

## References

  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Atri et al. (2020) Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2020, MNRAS, 493, L81, doi: 10.1093/mnrasl/slaa010
  * 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
  * Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65, doi: 10.1086/181711
  * Barnhart & Eichinger (2011) Barnhart, B. L., & Eichinger, W. E. 2011, Sol. Phys., 269, 439, doi: 10.1007/s11207-010-9701-6
  * 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
  * Blandford & Königl (1979) Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34, doi: 10.1086/157262
  * Bollimpalli et al. (2023) Bollimpalli, D. A., Fragile, P. C., & Kluźniak, W. 2023, MNRAS, 520, L79, doi: 10.1093/mnrasl/slac155
  * 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
  * Carson (1922) Carson, J. 1922, Proceedings of the Institute of Radio Engineers, 10, 57, doi: 10.1109/JRPROC.1922.219793
  * 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
  * 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
  * Daubechies et al. (2011) Daubechies, I., Lu, J., & Wu, H.-T. 2011, Applied and Computational Harmonic Analysis, 30, 243, doi: https://doi.org/10.1016/j.acha.2010.08.002
  * Dauser et al. (2016) Dauser, T., García, J., Walton, D. J., et al. 2016, A&A, 590, A76, doi: 10.1051/0004-6361/201628135
  * Dauser et al. (2013) Dauser, T., Garcia, J., Wilms, J., et al. 2013, MNRAS, 430, 1694, doi: 10.1093/mnras/sts710
  * Dauser et al. (2010) Dauser, T., Wilms, J., Reynolds, C. S., & Brenneman, L. W. 2010, MNRAS, 409, 1534, doi: 10.1111/j.1365-2966.2010.17393.x
  * 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
  * Denisenko (2018) Denisenko, D. 2018, The Astronomer’s Telegram, 11400, 1
  * 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
  * Esin et al. (1997) Esin, A. A., McClintock, J. E., & Narayan, R. 1997, ApJ, 489, 865, doi: 10.1086/304829
  * Fabian et al. (1989) Fabian, A. C., Rees, M. J., Stella, L., & White, N. E. 1989, MNRAS, 238, 729, doi: 10.1093/mnras/238.3.729
  * Fender (2001) Fender, R. P. 2001, MNRAS, 322, 31, doi: 10.1046/j.1365-8711.2001.04080.x
  * Fender et al. (2004) Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105, doi: 10.1111/j.1365-2966.2004.08384.x
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foster (1996) Foster, G. 1996, AJ, 112, 1709, doi: 10.1086/118137
  * Gao et al. (2023) Gao, C., Yan, Z., & Yu, W. 2023, MNRAS, 520, 5544, doi: 10.1093/mnras/stad434
  * García et al. (2014) García, J., Dauser, T., Lohfink, A., et al. 2014, ApJ, 782, 76, doi: 10.1088/0004-637X/782/2/76
  * García et al. (2016) García, J. A., Fabian, A. C., Kallman, T. R., et al. 2016, MNRAS, 462, 751, doi: 10.1093/mnras/stw1696
  * García et al. (2015) García, J. A., Steiner, J. F., McClintock, J. E., et al. 2015, ApJ, 813, 84, doi: 10.1088/0004-637X/813/2/84
  * Gilles (2013) Gilles, J. 2013, IEEE Transactions on Signal Processing, 61, 3999, doi: 10.1109/TSP.2013.2265222
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65, doi: 10.2140/camcos.2010.5.65
  * Guan et al. (2021) Guan, J., Tao, L., Qu, J. L., et al. 2021, MNRAS, 504, 2168, doi: 10.1093/mnras/stab945
  * 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
  * 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
  * Hu et al. (2022) Hu, C.-P., Lin, L. C.-C., Pan, K.-C., et al. 2022, ApJ, 935, 127, doi: 10.3847/1538-4357/ac8165
  * 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
  * Ingram & Done (2012) Ingram, A., & Done, C. 2012, MNRAS, 427, 934, doi: 10.1111/j.1365-2966.2012.21907.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 et al. (2015) Ingram, A., Maccarone, T. J., Poutanen, J., & Krawczynski, H. 2015, ApJ, 807, 53, doi: 10.1088/0004-637X/807/1/53
  * 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
  * Kalamkar et al. (2016) Kalamkar, M., Casella, P., Uttley, P., et al. 2016, MNRAS, 460, 3284, doi: 10.1093/mnras/stw1211
  * Kara et al. (2019) Kara, E., Steiner, J. F., Fabian, A. C., et al. 2019, Nature, 565, 198, doi: 10.1038/s41586-018-0803-x
  * Karpouzas et al. (2020) Karpouzas, K., Méndez, M., Ribeiro, E. M., et al. 2020, MNRAS, 492, 1399, doi: 10.1093/mnras/stz3502
  * Kato (1990) Kato, S. 1990, PASJ, 42, 99
  * Kawamuro et al. (2018) Kawamuro, T., Negoro, H., Yoneyama, T., et al. 2018, The Astronomer’s Telegram, 11399, 1
  * Koljonen et al. (2023) Koljonen, K. I. I., Long, K. S., Matthews, J. H., & Knigge, C. 2023, MNRAS, 521, 4190, doi: 10.1093/mnras/stad809
  * 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
  * Li et al. (2023) Li, X.-P., Yang, H.-Y., Cai, Y., et al. 2023, ApJ, 943, 157, doi: 10.3847/1538-4357/acae8c
  * 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
  * Lightman & Rybicki (1980) Lightman, A. P., & Rybicki, G. B. 1980, ApJ, 236, 928, doi: 10.1086/157820
  * 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., Ji, L., Bambi, C., et al. 2021, ApJ, 909, 63, doi: 10.3847/1538-4357/abdf65
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * 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
  * Markoff et al. (2005) Markoff, S., Nowak, M. A., & Wilms, J. 2005, ApJ, 635, 1203, doi: 10.1086/497628
  * 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
  * 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
  * Nathan et al. (2022) Nathan, E., Ingram, A., Homan, J., et al. 2022, MNRAS, 511, 255, doi: 10.1093/mnras/stab3803
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Shidatsu et al. (2019) Shidatsu, M., Nakahira, S., Murata, K. L., et al. 2019, ApJ, 874, 183, doi: 10.3847/1538-4357/ab09ff
  * Shui et al. (2023) Shui, Q. C., Zhang, S., Chen, Y. P., et al. 2023, 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 & Vietri (1998) Stella, L., & Vietri, M. 1998, ApJ, 492, L59, doi: 10.1086/311075
  * 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
  * 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
  * Sunyaev & Titarchuk (1980) Sunyaev, R. A., & Titarchuk, L. G. 1980, A&A, 86, 121
  * Tagger & Pellat (1999) Tagger, M., & Pellat, R. 1999, A&A, 349, 1003
  * Thomas et al. (2022) Thomas, J. K., Buckley, D. A. H., Charles, P. A., et al. 2022, MNRAS, 513, L35, doi: 10.1093/mnrasl/slab132
  * Titarchuk (1994) Titarchuk, L. 1994, ApJ, 434, 570, doi: 10.1086/174760
  * Torres et al. (2020) Torres, M. A. P., Casares, J., Jiménez-Ibarra, F., et al. 2020, ApJ, 893, L37, doi: 10.3847/2041-8213/ab863a
  * Tucker et al. (2018) Tucker, M. A., Shappee, B. J., Holoien, T. W. S., et al. 2018, ApJ, 867, L9, doi: 10.3847/2041-8213/aae88a
  * Uttley et al. (2018) Uttley, P., Gendreau, K., Markwardt, C., et al. 2018, The Astronomer’s Telegram, 11423, 1
  * 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
  * Veledina et al. (2013) Veledina, A., Poutanen, J., & Ingram, A. 2013, ApJ, 778, 165, doi: 10.1088/0004-637X/778/2/165
  * Wagoner (1999) Wagoner, R. V. 1999, Phys. Rep., 311, 259, doi: 10.1016/S0370-1573(98)00104-5
  * Wang et al. (2022) Wang, J., Kara, E., Lucchini, M., et al. 2022, ApJ, 930, 18, doi: 10.3847/1538-4357/ac6262
  * Wijnands et al. (1999) Wijnands, R., Homan, J., & van der Klis, M. 1999, ApJ, 526, L33, doi: 10.1086/312365
  * Yang et al. (2022) Yang, Z.-X., Zhang, L., Bu, Q.-C., et al. 2022, ApJ, 932, 7, doi: 10.3847/1538-4357/ac63af
  * You et al. (2018) You, B., Bursa, M., & Życki, P. T. 2018, ApJ, 858, 82, doi: 10.3847/1538-4357/aabd33
  * You et al. (2020) You, B., Życki, P. T., Ingram, A., Bursa, M., & Wang, W. 2020, ApJ, 897, 27, doi: 10.3847/1538-4357/ab9838
  * You et al. (2021) You, B., Tuo, Y., Li, C., et al. 2021, Nature Communications, 12, 1025, doi: 10.1038/s41467-021-21169-5
  * Yu et al. (2023) Yu, W., Bu, Q.-C., Yang, Z.-X., et al. 2023, arXiv e-prints, arXiv:2305.12317, doi: 10.48550/arXiv.2305.12317
  * Zdziarski et al. (1996) Zdziarski, A. A., Johnson, W. N., & Magdziarz, P. 1996, MNRAS, 283, 193, doi: 10.1093/mnras/283.1.193
  * Zdziarski et al. (1999) Zdziarski, A. A., Lubiński, P., & Smith, D. A. 1999, MNRAS, 303, L11, doi: 10.1046/j.1365-8711.1999.02343.x
  * 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. (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
  * Zhao et al. (2021) Zhao, X., Gou, L., Dong, Y., et al. 2021, ApJ, 916, 108, doi: 10.3847/1538-4357/ac07a9
  * Zhou et al. (2022) Zhou, D.-K., Zhang, S.-N., Song, L.-M., et al. 2022, MNRAS, 515, 1914, doi: 10.1093/mnras/stac1789
  * Życki et al. (1999) Życki, P. T., Done, C., & Smith, D. A. 1999, MNRAS, 309, 561, doi: 10.1046/j.1365-8711.1999.02885.x