A Phase-resolved View of the Low-frequency Quasiperiodic Oscillations from the Black Hole Binary MAXI J1820+070

Although low-frequency quasiperiodic oscillations (LFQPOs) are commonly detected in the X-ray light curves of accreting black hole X-ray binaries, 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 quasiperiodic oscillation (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.


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 (∼ 100 keV) (Sunyaev & Titarchuk 1980;Titarchuk 1994;Zdziarski et al. 1996Zdziarski et al. , 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α lines at ∼ 6.4 keV and a Compton hump in the ∼20-30 keV energy range, etc. (Lightman & Rybicki 1980;Fabian et al. 1989;Dauser et al. 2010;García et al. 2014, 2015, andreferences 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 ∼ 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. 2021Ma et al. , 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. 2016Ingram et al. , 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 Xray transient by MAXI/GCS on 2018 March 11 (Denisenko 2018;Kawamuro et al. 2018;Tucker et al. 2018), located at R.A.=18 h 20 m 21 s .9 and decl.=+07• 11 ′ 07 ′′ .3(Kawamuro et al. 2018).So far, it is one of the brightest X-ray transients, having a soft X-ray flux of ∼ 4 Crab in 2-4 keV (Shidatsu et al. 2019) with a low column density of 1.5 × 10 21 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 ± 0.33 kpc.The VLBI observation provides an inclination angle of 63 ± 3 • and a mass of 9.2 ± 1.3M ⊙ (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 ⊙ .Guan et al. (2021) and Zhao et al. (2021) estimated the BH spin as 0.2 +0.2 −0.3 and 0.14 ± 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.

OBSERVATIONS AND DATA REDUCTION
In this study, we analyze data collected by the Hard X-ray Modulation Telescope (Insight-HXMT) (Zhang et al. 2014(Zhang et al. , 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 ∼ 0.04 to ∼ 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,209,MJD 58,225,MJD 58,239,and MJD 58,253,respectively.From Epochs 1 to 4, the count rate of LE evolves from ∼ 1000 to ∼ 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.06 1 and the standard Insight-HXMT Data Reduction Guide v2.05 2 , and apply the following series of criteria recommended by the Insight-HXMT team: (1) the elevation angle (ELV) is greater than 10 • ; (2) the geometric cutoff rigidity (COR) is greater than 8 GeV; (3) the pointing position offset is less than 0.04 • ; 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).

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: here the phase of the k-th mode, denoted as ϕ k (t), is a non-decreasing function, indicating that the instantaneous frequency, ω k (t) = ϕ ′ k (t) ≥ 0. Additionally, the envelope A k (t) is non-negative (A k (t) ≥ 0), while both the instantaneous frequency ω k (t) and envelope A k (t) exhibit much slower variations compared to the phase ϕ 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 ω k .For a mode, u k (t), with a mean frequency, w k , its practical bandwidth increases both with the maximum deviation ∆f of the instantaneous frequency from its center, and with the rate of this excursion, f 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 AM , broadens the spectrum even further.Therefore, the total practical bandwidth of an IMF can be estimated as: 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 gradient3 : where BW k indicates the bandwidth of the mode, δ 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 {ω k } = {ω 1 , ..., ω K }.The objective is to minimize the sum of bandwidths of the modes by solving the following optimization problem: 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 α serves as a weighting factor that balances the compactness and reconstruction error.The VMD algorithm iteratively updates {u k } and {ω k } until convergence.In this process, ω k is directly calculated in the Fourier domain as: where ûk (ω) 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 α is crucial in the VMD algorithm and determines its performance.Typically, higher values of K are associated with higher values of α, indicating that as the number of modes increases, the bandwidth of each mode decreases.In this study, our ultimate goal is to perform the phaseresolved 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 α 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 α 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.24 (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 ∼0.3 Hz.In this example, we set K = 5 and α = 4000 for the VMD algorithm.A ∼ 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.

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 In this example, we set K = 5 and α = 4000 for the VMD algorithm, hence obtain five IMFs.The IMF C2 corresponds to the QPO signal at ∼ 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.where p.v. is the Cauchy principal value.Using the Hilbert transform, one can obtain the corresponding analytic signal, which is defined as where j is the imaginary unit.The time-dependent amplitude, A(t), and phase, ϕ(t), can be obtained from respectively.Therefore, the instantaneous frequency, ν(t), can be defined as 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 ∼0.3 Hz oscillation is relatively weak in time intervals of ∼ 18-25 s and ∼ 65-80 s.However, it becomes stronger during time intervals of ∼ 40-60 s and ∼ 80-90 s, with frequencies of ∼ 0.27 and ∼ 0.35 Hz, respectively (see Figure 4).Once the instantaneous phases, ϕ(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 ∼0.7 cycles, followed by a rapid drop of ∼0.3 cycles.This non-sinusoidal behavior is associated with the presence of the second harmonic in the PDS (see Figure 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×T babs×(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 T babs accounts for interstellar absorption, with the equivalent hydrogen column, N H , fixed at 0.15 × 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 in ) and the inner radius (R in ) of the accretion disk.The latter one is linked to the normalization of the diskbb model follows as N disk = (R in /D 10 ) 2 cos i, where N disk is the normalization, R 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 ⊙ (Torres et al. 2020), the distance and inclination as 2.96 kpc and 63 • , respectively (Atri et al. 2020), we can obtain N disk ≈ 816(R in /R g ) 2 , where R g = GM/c 2 is the gravitational radius.
For the relativistic reflection model relxillCp, we assume the canonical power-law emissivity profile, ϵ = 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 • , 0.14, respectively (Atri et al. 2020;Zhao et al. 2021).The parameter R 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, 1000R g .Since R in can also be a free parameter in relxillCp, it is self-consistent to link it to N disk in diskbb as N disk = 816R 2 in (see also Gao et al. 2023), where R in is in units of R 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 e ) is set to the default value of 10 15 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 ξ = 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 Fe and the ionization parameter log ξ to the same values obtained from the phase-averaged spectral fittings.After fitting the energy spectra with the model, we obtain a satisfactory reduced χ 2 (∼ 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 displays the phase dependence of the seven free parameters, including the temperature of the inner radius of the disk (T in ), inner radius (R in ), spectral index (Γ), electron temperature (kT e ), reflection fraction (R 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.From the phase-resolved spectral analysis, we observe no significant modulations in T in and R 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 f and kT e in Epochs 3 and 4, but not 1 and 2. For the modulation of R f in Epochs 3 and 4, it seems to lead the flux by ∼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.

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 phaserevolved 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.

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 (ψ/π) 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 ψ/π 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 ψ/π distributions, but the HHT method produces a more concentrated distribution.The mean values of ψ/π from the Fourier and HHT methods are found to be 0.191 ± 0.041 and 0.229 ± 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, ψ/π, 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.
de Ruiter et al. ( 2019) measured the phase difference for many observations of QPOs in the RXTE archive.In most high-inclination systems, ψ/π between harmonics is distributed around ∼ 0 to ∼ 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 ψ may slightly alter the waveform.Considering that the inclination of MAXI J1820+070 is ∼ 63 • , 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 ψ/π ∼ 0.2 for QPOs at ∼ 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 ψ in MAXI J1820+070 with that in other sources, it is necessary to have more QPO samples at higher QPO frequencies (∼ 1 − 5 Hz) which were not detected by Insight-HXMT.

Modulation of the Spectrum
As shown in Figure 6, the spectral index (Γ) and the normalization of relxillCp (Norm 1 ) 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. 2021Ma et al. , 2023)).

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(Ingram et al. , 2015;;You et al. 2018You et al. , 2020;;Shui et al. 2023).
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 e ) and the optical depth (τ ).Therefore, the modulations of Γ could be attributed to the modulations of kT e and τ .For a given kT e and spectral index (Γ), we can estimate the corona's optical depth by where m 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 Γ and kT 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: where τ 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.θ is the included angle between the line of sight and the corona's normal, and it varies over the precession period.The cosine of θ can be expressed using the following equation: where i is the inclination angle between the line of sight of the observer and the black hole spin axis, Φ is the azimuth angle of the observer, γ is the included angle between the corona's normal and the black hole spin axis and ω is the precession phase angle, which varies from 0 to 2π during the precession period.We set i = 63 • , Φ = 180 • (resulting in the peak flux occurring at 0.5 QPO cycles) and γ = 5 • 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 τ 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 τ .
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 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 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. 2018You et al. , 2020;;Nathan et al. 2022), subsequently leading to corresponding modulations in the derived reflection fraction (see You et al. 2020).

Jet Precession Model
Another geometric model worth considering is the jet precession model proposed by Ma et al. (2021Ma et al. ( , 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 Γ 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 ∼ 100 keV with respect to that at ∼ 2 keV) are ∼0.5 s at the QPO frequency of ∼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 where F obs is the observed spectrum, F int is the intrinsic spectrum of the jet, Γ is the spectral index of the jet emission and D is the relativistic boosting factor given by where β = v/c is the jet speed in units of the speed of light, and θ 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 where F p and F t are the spectra at the peak and trough phases, respectively, and θ p and θ t are the included angles at the peak and trough phases, respectively.In the case of Φ = 180 • , θ t = θ ω=0 and θ p = θ ω=π , respectively.We define the ratio between the trough spectrum and peak spectrum as R = F t /F 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 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 where E is the energy of emitted photons.Here, we set Γ as the value obtained from the spectral fitting of the phaseaveraged spectra.The bottom panels of Figure 10 displays the energy dependence of β.The obtained anti-correlation between β 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 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 (γ).In fact, the coupling between the jet and inner disk (see Ma et al. 2023) could entail a decrease in γ 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 γ.Given that the QPO RMS is influenced by both the jet speed and γ, the combined impact of jet acceleration and the γ 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 β-E correlation within each epoch remains unaffected by the assumption of a constant γ.
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 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 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 χ 2 (∼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 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 (β) across the four epochs, where values of β are calculated by Equation 16.The gray scatters represent the calculation from the original spectra, while colored scatters represent rebinned values.
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 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 magnetohydrodynamic (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 Norm 2 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 ∼ 20 s, the corresponding light-crossing distance is estimated to be ∼ 6 × 10 6 km (equivalent to ∼ 4.8 × 10 5 R g , assuming M BH = 8.5M ⊙ ).Based on the well determined orbital parameters of the system, the outer disc radius is constrained to be between ∼ 6.3 × 10 4 and ∼ 2.5 × 10 5 R 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 Norm 2 in Epochs 1 and 2 is not likely from reflection from the distant disk.

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).

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). .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.

Figure 1 .
Figure1.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.

Figure 2 .
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 α = 4000 for the VMD algorithm, hence obtain five IMFs.The IMF C2 corresponds to the QPO signal at ∼ 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 4 .
Figure 4. Left: Hilbert spectrum of ∼ 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.

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

Figure 6 .
Figure6.Phase dependence of the seven free spectral parameters, including the temperature of the inner radius of the disk (Tin), inner radius (Rin), spectral index (Γ), electron temperature (kTe), reflection fraction (R 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 7 .Figure 8 .
Figure7.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 π. Right: Reconstructed QPO waveform derived the HHT phase-folding (black) and Fourier (red) methods.

Figure 9 .
Figure9.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 11 .
Figure11.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.
Figure12.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.

Table 1 .
Log of Insight-HXMT Observations Used in This Work