Recovery of High-energy Low-frequency Quasiperiodic Oscillations from Black Hole X-Ray Binary MAXI J1535–571 with a Hilbert–Huang Transform Method

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

Following a prolonged period of quiescence at lowflux levels, black hole X-ray binaries (BHXRBs) can experience outbursts, with substantial changes in their X-ray emission properties (Remillard & McClintock 2006;Done et al. 2007).During a complete outburst, the source typically exhibits four distinct canonical states: the low/hard state (LHS), hard intermediate state (HIMS), soft intermediate state (SIMS) and high soft state (HSS).A source exhibits distinct X-ray spectral and variability properties in different states (Belloni et al. 2005;Homan & Belloni 2005).Low frequency quasi-periodic oscillations (LFQPOs) are observed in most transient BHXRBs, typically ranging from 0.1 to 30 Hz (van der Klis 1989).These oscillations are characterized by a narrow peak with finite width in the power density spectra (PDS) and are typically classified into three types (A, B, and C) based on their centroid frequency, Q factor (ratio of frequency to the full width at half maximum), and root-mean-square (rms) amplitude (Wijnands et al. 1999;Casella et al. 2005;Remillard & McClintock 2006).
LFQPOs are thought to originate from the innermost part of the accretion flow, where the most X-ray emissions and fast variability are produced.Fast X-ray variability is observed across a wide range of timescales (Belloni & Stella 2014).In addition to the QPO component, broadband continuum components also contribute to the PDS, which represent fast aperiodic variability.Several models have been suggested to account for this noise component, such as the shot noise model, coronal flare model, and fluctuation propagation model (Terrell 1972;Nolan et al. 1981;Scargle et al. 1993;Mineshige et al. 1994;Ingram & Done 2011).For LFQPOs, various models have been also proposed, based on either the geometric effect or the intrinsic variability of the accretion flow.Intrinsic models include the trapped corrugation modes (Kato 1990;Wagoner 1999), the Accretion-ejection instability model (AEI, Tagger & Pellat 1999), and the Two-Component Advection Flow model (Molteni et al. 1996), etc.Most geometric models, on the other hand, are associated with the relativistic Lense-Thirring (L-T) precession effect, such as the relativistic precession model (RPM; Stella et al. 1999), L-T precession of the entire hot accretion flow (Ingram et al. 2009) and precession of the jet (Stevens & Uttley 2016;de Ruiter et al. 2019;Ma et al. 2021Ma et al. , 2023)).
Componization is widely considered to be the radiative mechanism behind type-C QPOs (e.g., Méndez et al. 2022;Bellavita et al. 2022), because observational studies show that the factional rms of type-C QPOs generally increases with photon energy (Zhang et al. 2017;Huang et al. 2018;Kong et al. 2020;Zhang et al. 2020a) and no evident disc-like component is detected in their rms spectra (Sobolewska & Życki 2006;Axelsson et al. 2014;Axelsson & Done 2016;Shui et al. 2023a).Furthermore, Ma et al. (2021) reported the discovery of tpye-C QPOs above 200 keV, further indicating a strong association with Comptonized emission.Additionally, inclination-dependent amplitudes and time lags of type-C QPOs (see Motta et al. 2015;van den Eijnden et al. 2017), and QPO phase-dependent reflection component revealed with phase-resolved spectroscopy (see Ingram & van der Klis 2015;Ingram et al. 2016Ingram et al. , 2017;;Nathan et al. 2022;Shui et al. 2023a), may suggest type-C QPOs originate from some geometrical effects.For recent reviews of observations and theories of LFQPOs, the reader is referred to Ingram & Motta (2019).
The study of LFQPOs in high energy bands (> 30 keV) is essential to our understanding of their origin, since the radiative mechanism and geometric properties of LFQPOs could be further constrained (Huang et al. 2018;Kong et al. 2020;Ma et al. 2021;Liu et al. 2021).Thanks to the wide energy range of Insigt-HXMT, we can study high-energy LFQPOs above 100 keV in some BHXRBs in details (see e.g.Ma et al. 2021).Furthermore, γ-ray QPOs have been detected by Fermi-LAT in numerous bright blazars (see e.g.Sandrinelli et al. 2014;Ackermann et al. 2015;Banerjee et al. 2023).Consequently, certain BHXRBs may have the potential to exhibit γ-ray QPOs, which could involve different radiative mechanisms from those observed in the X-ray bands, and warrant further detailed investigation in the future.In this Letter, we propose a novel technique based on the Hilbert-Huang transform (HHT) for recovering highenergy QPO signals in MAXI J1535-571, which is one of the brightest BHXRBs observed by Insight-HXMT.The HHT is a powerful tool for analyzing phenomena with non-stationary periodicity proposed by Huang et al. (1998) and has been successfully applied in astronomical researches (Camp et al. 2007;Hu et al. 2011Hu et al. , 2014;;Su et al. 2015;Hsieh & Chou 2020;Yu et al. 2023a;Shui et al. 2023a).The HHT allows us to track changes in the frequency of QPOs on short time-scales and conduct their phase-resolved spectral analysis.
In this Letter, we propose a HHT-based technique to reconstruct the high-energy QPO waveforms.In Section 2, we outline the Insight-HXMT observations and data reductions.Then, in Section 3, we present a brief introduction to the HHT method and its application in recovering high-energy QPO signals, and then present the corresponding results derived from MAXI J1535-571.Finally, we discuss and summarize these results in Section 4 and Section 5, respectively.
We focus on high-energy signals of the type-C QPO detected by Insight-HXMT, so one of the bright systems MAXI J1535-571 with numerous type-C QPO observations (see Huang et al. 2018;Kong et al. 2020;Zhang et al. 2022;Yu et al. 2023b) is investigated in this study.We extract the event data using Insight-HXMT Data Analysis Software v2.05 1 and the current calibration model v2.06 2 , following the criteria recommended by the Insight-HXMT team: (1) elevation angle (ELV) greater than 10 • ; (2) geometric cutoff rigidity (COR) greater than 8 GeV; (3) pointing position offset less than 0.04 • ; and (4) good time intervals (GTIs) at least 300 s away from the South Atlantic Anomaly (SAA).Backgrounds are derived from blind detectors using the LEBKGMAP, MEBKGMAP, and HEBKGMAP tools, version 2.0.9, based on the standard Insight-HXMT background models (Liao et al. 2020a,b;Guo et al. 2020).All photons are corrected to the barycenter of the solar system for their arrival times using the HXMTDAS tool hxbary.In Figure 1, we present the long-term net light curve for MAXI J1535-571 obtained from HE in the energy range of 30-120 keV.To ensure robust count statistics for HE, we select QPO observations with HE count rates in the energy range of 30-120 keV that closely approximate the highest level across the outburst (as indicated by the shaded area in Figure 1).We subsequently generate net light curves with a time resolution of 0.05 s for these selected observations, which are then subjected to  the detailed timing analysis.Detailed information regarding the aforementioned observations is available in Table 1.

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

Hilbert-Huang Transform Analysis
The HHT, proposed by Huang et al. (1998) as an adaptive data analysis method, offers a powerful tool for analyzing signals with non-stationary periodicity.It consists of two main components: mode decomposition, which aims to decomposing a signal into several intrinsic mode functions (IMFs), and Hilbert spectral analysis (HSA), which allows for the extraction of frequency and phase functions of the desired IMFs, such as the QPO component (e.g.Hu et al. 2014;Hsieh & Chou 2020).Following Shui et al. (2023a), we use the variational mode decomposition (VMD, Dragomiretskiy & Zosso 2014) method to extract the IMF correspond to the QPO component.Compared to the traditional method empirical mode decomposition (EMD, Huang et al. 1998;Huang & Wu 2008) technique, the VMD method theoretically mitigates the influence of mode mixing issue by decomposing the signal into a sum of IMFs with analytically calculated limited center frequency and bandwidth.Mode mixing is one of the major drawbacks of EMD, which is defined as a single IMF either consisting of signals of widely disparate scales or a signal of a similar scale residing in different IMF components (Huang & Wu 2008).In the VMD algorithm, the set of IMFs of an input signal, f (t), is initialized as {u k } = {u 1 , ..., u K } along with their corresponding center frequencies {ω k } = {ω 1 , ..., ω K }.The objective of the VMD is to minimize the sum of bandwidths of IMFs by solving the following optimization problem 3 : where K is the total number of modes, BW k is the bandwidth of the kth mode.The first component of the objective function enhances mode compactness, whereas the second component quantifies the reconstruction error, which ought to be minimized.The parameter α acts as a balancing weight between compactness and the reconstruction error.In Appendix A, we offer a concise overview of the VMD algorithm.For the complete algorithmic description of the VMD, we recommend referring to Dragomiretskiy & Zosso (2014).There are two parameters should be initially set in the VMD algorithm, which is the total number of modes K and the weighting factor α. In the VMD algorithm, selecting the parameters K and α is vital for its efficacy.Generally, an increase in K correlates with a rise in α, suggesting that the bandwidth narrows as the number of modes grows.Following the procedure described in Shui et al. (2023a), we initially assign α a typical value of 2000 and opt for 2 modes (K = 2).We then generate the power spectrum for each mode and compare it with the power spectrum of the original light curve.Should no mode display a power spectrum mirroring the QPO distinctive . From another perspective, the inner product can be expressed as ⟨p(x), q(x)⟩ = p * (x)q(x)dt, which implies that ||h(x pattern, we increment K.This iterative strategy persist until the fundamental and harmonic components of the QPO are discernibly isolated in individual modes.Finally, we adjust α to align the IMF PDS profile with the QPO peak in the power spectrum of the original light curve.The selection of K and α involves elements of subjectivity, which is one of the major concerns of the VMD algorithm (Nazari & Sakhaei 2018).This issue requires future investigations, which are beyond the scope of this study.Nevertheless, we find that making moderate adjustments to the values of K and α have a minimal impact on the subsequent Hilbert analysis when ensuring that the PDS of the desired IMF aligns with the QPO peak in the observed PDS (see Figure 2).In addition, the current VMD algorithm addresses certain limitations of the original EMD algorithm, and show attractive performance compared to the existing mode decomposition methods.The code we used to perform VMD is from vmdpy v0.24 (Carvalho et al. 2020), which is an open-source Python package.It is worth noting that an appropriate time resolution is required for the HHT analysis, as too dense sampling results in poor SNR in each bin, while too sparse sampling does not sufficiently decompose the QPO component.In our case, where the QPO frequency is ∼2-3 Hz and the net count rate is ∼ 600 cts s −1 , a time resolution of 0.05s is appropriate.This ensures that the time internal between the samplings is much shorter than the variability timescale, and each time bin has a signal-to-noise ratio (SNR) ≳ 5σ. Figure 2a shows an example of a 14-s-long light curve with a 0.05s time resolution from HE observations of MAXI J1535-571 (obsID P011453500201) and the corresponding IMFs.In this representative observation, the observed light curve exhibits an LFQPO around 2 Hz.
For the VMD algorithm, we selected parameters with K = 5 and α = 300.We successfully identify a ∼2 Hz oscillation and its harmonic component as the second and third IMFs, respectively.To compare results of the VMD algorithm with the Fourier analysis, Figure 2b plots the power spectra for the original HE light curve alongside that of the second and third IMFs.These power spectra are generated in the frequency range of 0.04-10 Hz.A simulation also demonstrates the robustness of the VMD algorithm that it shows good performance in the presence of noise in the analyzed observations (see Appendix B for details).
For each IMF obtained from the mode decomposition, the Hilbert transform can calculates its physically meaningful instantaneous phase, amplitude and frequency functions.For a time series f (t), its Hilbert transform is defined as 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 only apply the Hilbert transform to the IMF corresponding to the QPO component (e.g. the second IMF plotted in Figure 2a).From the top to bottom panels of Figure 2c, we illustrate the QPO IMF and corresponding phase function, instantaneous frequency and amplitude of the QPO component, respectively.As one can see, the instantaneous frequency of QPO fluctuates between ∼1 and ∼3 Hz over time.

Recovery of QPO Waveforms above 120 keV
We perform the HHT analysis for each selected observations from HE/NaI in the energy range of 30-120 keV.In such an energy band, HE shows a good performance to provide high SNR light curves.These high SNR light curves are suitable to be performed with the HHT analysis.Once the instantaneous phases, ϕ(t), of the QPOs have been computed using the Hilbert transform as plotted in Figure 2c, the QPO waveforms can be constructed by phase-folding the light curves (also see Hu et al. 2014;Hsieh & Chou 2020;Shui et al. 2023a).We plot the reconstructed 30-120 keV QPO waveform in the top panel of Figure 3.To enhance the SNR, waveforms from all selected observations have been collectively combined, and the final combined waveform are normalized by the mean value of the count rate.It is evident that the QPO waveform exhibits an non-sinusoidal nature, which is characterized by a fast rise of ∼ 0.4 cycles, followed by a gradual drop of ∼ 0.6 cycles.Utilizing the well-defined QPO phases derived from the HHT analysis, we extend the analysis to perform phase-folding on high-energy light curves above 120 keV.It should be noted that, due to the low SNR in the high-energy band, conducting the HHT analysis there is not feasible.Consequently, the phase-folding of the high-energy light curves is carried out using the instantaneous phase obtained from the HHT analysis in the low-energy band, assuming that the high energy QPOs have the same frequency as that in the low energy band.From the second to the bottom panels of Figure 3, the phase-folded light curves in various energy bands above 120 keV are plotted.We use two methods to determine the significance of the high energy QPO profiles, one is the standard χ 2 -test method, another is the crosscorrelation technique described in Hou et al. (2022).It is evident that the QPO profile remains distinctly identifiable in the energy bands between 120 and 200 keV.However, signals dissipate in energy ranges exceeding 200 keV.The significance of the QPO profile in the 150-170 keV energy band derived from the χ 2 -test method is 5.2σ, while the cross-correlation technique gives a 8.1σ confidence level.As for the QPO profile in 170-200 keV, the cross-correlation technique gives a 6.2σ confidence level (see Appendix C for details).However, the χ 2 -test method gives only 2.1σ for the 170-200 keV profile.This indicates that the cross-correlation method is more advantageous since it takes the lower-energy QPO profile, which has a very high SNR, as a template in searching for the higher-energy pulsed signals (see also Hou et al. 2022).We also plot the PDS for high-energy light curves (100-200 keV) to provide a comparison between the HHT phase-folding and Fourier methods.Notably, no significant QPO peaks are observed in the long-term combined PDS above 100 keV, as shown in Figure 4.In a simulation study presented in Appendix B, we find that a small part of the noise component is decomposed into the QPO IMF, resulting from mode mixing.How-0 .9 0 1 .0 5 1 .2 0 0 .9 8 0 1 .0 0 8 1 .0 3 6 0 .9 8 1 .0 0 1 .0 2 0 .9 9 0 1 .0 1 2 0 .9 9 1 6 0 .9 9 8 3 1 .0 0 5 0 0 .9 9 6 4 1 .0 0 1 7 1 .0 0 7 0 0 0 .5 1 1 .5 2 0 .9 9 5 4 ever, the simulation results show that if the QPO signals do not exist in the 170-200 keV energy band, the significance of the modulation from the noise is found to be < 2σ, which is inconsistent with the observational results (∼ 6.2σ).Conversely, when QPO signals are added to the high-energy light curves, the significance of the modulation in the phase-folded light curve is found to be ∼ 6 − 8σ, consistent with the observational results.This provides supports for the idea that the observed modulations of the phase-folding profile in the energy ranges of 170-200 keV could indeed originate from the QPO component, despite the presence of some mode mixing issues in the VMD.(see Appendix B for details).
F r e q u e n c y ( H z )

Properties of High Energy QPOs
It is crucial to highlight that the 170-200 keV QPO signal, as reported, marks the highest-energy QPO detection for this source to date.To further study the properties of the high-energy (>120 keV) QPO profiles, we analyze the phase-folded waveforms using harmonically related cosine functions for fitting.From Figure 3, it is evident that the QPO waveforms exhibit a nonsinusoidal nature, particularly below 130 keV, which is consistent with the emergence of the second harmonic in the PDS (see Figure 2b).We find that two cosine functions provide a good description of the data.Furthermore, we subtract the background contribution from the fractional amplitude, resulting in the final fitting function of where A 1 and A 2 represent the amplitudes, ϕ 1 and ϕ 2 the phases of the cosine functions, and S and B denote the mean count rates of the source and background, respectively.We perform phase-folding and subsequent model fitting in 8 distinct energy bands of 30-40, 40-50, 50-70, 70-100, 100-130, 130-150, 150-170, and 170-200 keV.The analysis excludes waveforms above 200 keV due to the insignificance of the QPO signal in these bands.Notably, the parameters of the second harmonic component (A 2 and ϕ 2 ) are effectively constrained only at energies below 130 keV.This suggests a notable attenuation of the second harmonic component above 130 keV, prompting us to employ one cosine function for fitting in these energy bands.Figure 5 presents the energy dependence of the best-fitting parameters.It is clear that ϕ 1 and ϕ 2 remains relatively stable from 30 to 130 keV, whereas ϕ 1 exhibits a decreasing trend with increasing energy beyond 130 keV, indicating soft lags across different energy bands.Regarding the amplitudes, while A 2 stays relatively constant, A 1 shows a slight increase from 30 to 130 keV, followed by a decrease from ∼ 0.3 to ∼ 0.1 above 130 keV.

DISCUSSION
In this study, we have performed a comprehensive analysis of type-C QPOs in the black hole binary MAXI J1535-571, utilizing data from the Insight-HXMT.Based on the Hilbert-Huang transform (HHT) method, we have proposed a novel method for detecting high-energy QPO signal.This allows us to successfully reconstruct the QPO waveform above 100 keV through phase-folding of the light curves, although no significant QPO signal has been discernible in the PDS beyond this threshold.Notably, by combining QPO waveforms derived from seven observations, we have achieved a highly confident detection of QPO profiles above 170 keV at a significance level of ∼ 6.2σ.This represents the highestenergy detection of QPOs in this source reported to date.Furthermore, this technique has the potential to extend the study of energy-dependent behaviors of QPO amplitude and phase lag up to 200 keV.
The general method for detecting QPO signals from X-ray light curves is the Fourier transform, which allows QPOs to be identified as a series of harmonic peaks in the PDS.However, in the high-energy bands (>100 keV), the limitations in count statistics present a challenge.Previous studies commonly combined multiple observations to enhance the SNR of QPO peaks in the PDS (see e.g.Ma et al. 2021;Yang et al. 2022).However, this approach assumes similar QPO frequencies across the combined observations, which is not typically the case during the HIMS of outburst rising phase of BHXRBs, where the QPO frequency tends to increase rapidly (see e.g.Huang et al. 2018;Shui et al. 2021Shui et al. , 2023b)).If the QPO frequencies among combined observations show significant differences, multiple fundamental QPO peaks would appear in the long-term combined PDS in the low-energy band (see Figure 4).In such cases, improving the SNR of higher energy QPO signals would be challenging, as the QPO signals from different observations are distributed over a wide frequency range.Although any method would face with the limitations in count statistics in the high-energy bands, the HHT phase-folding method proposed in this work effectively addresses the aforementioned issue encountered in the Fourier method when combining different observa-tions.In the phase-folding approach, the well-defined QPO phases can be obtained in the low-energy band.This indicates that, for any QPO frequency, a QPO cycle can be divided into given phase bins between 0 and 2π.Consequently, by stacking up counts within the same phase bin from multiple observations, one can achieve a higher SNR associated with the modulation of the count rate as function of the phase (i.e. the QPO waveform).Both our observational and simulation studies demonstrate that, provided the phase information of the QPO, the phase-folding approach performs better in determining the high-energy QPO signal than the Fourier method when the QPO frequency varies among different observations (see Appendix B).
Previous studies by Huang et al. (2018) and Kong et al. (2020) have presented the energy dependence of type-C QPO amplitude in MAXI J1535-571 across a wide energy band ranging from 1 to 100 keV.They observed an increase in QPO RMS amplitude with photon energy up to ∼20 keV, after which it plateaus across all observations.Utilizing our HHT-based technique, we extend this study to analyze QPO amplitude beyond 100 keV.As shown in Figure 5b, we note a marginal increase in amplitude from 30 to 60 keV, followed by a plateau between 60 and 120 keV, and a decreasing trend beyond 120 keV.The flat trend of amplitude between 30 and 120 keV agrees with findings reported by Huang et al. (2018) and Kong et al. (2020).
In broad-band spectral analyses of MAXI J1535-571, both Kong et al. (2020) and Rodi et al. (2022), using Insight-HXMT and INTEGRAL/SPI data respectively, reported findings indicative of a relatively lowtemperature corona.Kong et al. (2020) reported the cutoff energy of the power-law component within the 60-80 keV range, while Rodi et al. (2022) found the temperature of thermal electrons in the corona to be ∼ 20 keV during the HIMS.These observations suggest that emissions above ∼100 keV are unlikely to be coronadominated.Furthermore, using INTEGRAL/SPI data, Rodi et al. (2022) identified an additional high-energy component, which can be modeled by a cutoff powerlaw and is observed to have a flux approximately six times over the coronal flux in the 100-400 keV range (refer to Figure 6 in Rodi et al. (2022)).Consequently, it is probable that the QPO amplitude above 100 keV is predominantly influenced by this extra high-energy component.The different evolutionary trends between the coronal and extra power-law fluxes, as reported by Rodi et al. (2022), imply different origins for these two components, with the high-energy component likely being from a jet.This is further corroborated by the divergent energy-dependent properties of both ϕ 1 and A 1 in low- (30-130 keV) and high-energy (130-200 keV) bands, as shown in Figure 5. Particularly, the emergence of soft lags above 130 keV agrees with predictions made by the jet precession model (Ma et al. 2021(Ma et al. , 2023;;Shui et al. 2023a).
Using the phase function in 30-120 keV to fold the light curves above 120 keV appears to introduce some systematic issues, as spectral analyses have indicated different dominant components between the two energy bands.However, our results demonstrate that the highenergy QPO waveform can be effectively reconstructed using the 30-120 keV phase function, suggesting correlated QPO phases across these energy bands.We propose that a jet-corona coupling precession scenario could explain our findings, wherein the jet and corona are coupled and precess with the same period, leading to correlated QPO signals from the two locations.This picture is further supported by a recent general relativistic magneto-hydrodynamic simulation (Liska et al. 2018).
Given that typical speeds at the base of the jet can already be relativistic (McKinney 2006), a fundamental concept in the jet precession model is that emission along the jet axis is considerably stronger than off-axis emission due to the Doppler boosting effect (see Ma et al. 2021).As a result, when the jet precesses, the flux from it undergoes modulation.The velocity of the jet material significantly influences the QPO variability.In particular, higher jet speeds result in greater modulation amplitude, as there is a larger difference in emission strength between on-axis and off-axis situations.During the stage of the studied observations, it was reported to exist a compact and continuous jet in MAXI J1535-571 from the radio band (Russell et al. 2019).For a jet of this nature, the observed luminosity, S obs , is described by where S obs represents the observed luminosity, S int is the intrinsic luminosity of the jet, Γ is the spectral index of the power-law spectrum 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 angle between the jet axis and the observer's line of sight, calculated as cos θ = sin i cos Φ sin γ cos ω + sin i sin Φ sin γ sin ω where i is the inclination angle between the observer's line of sight and the black hole spin axis, Φ is the azimuth angle of the observer, γ is the angle between the jet axis and the black hole spin axis, and ω is the precession phase angle, ranging from 0 to 2π over the precession period.(see Ma et al. 2021, for details).Setting the inclination, i, to 57

SUMMARY AND CONCLUSION
A detailed analysis of type-C QPOs from the black hole binary MAXI J1535-571 with observations from Insight-HXMT is presented in this work.We propose a novel method based on the Hilbert-Huang transform, which enables a robust detection of type-C QPOs up to 200 keV, and consistently above 170 keV.This method also has the potential to extend the study of the energy dependence of QPO amplitude and phase lag up to 200 keV.Our findings reveal that the QPO fundamental exhibits soft lags and a decrease in amplitude above ∼ 100 keV.Given the reported detection of a hard tail in broad spectroscopy, we proposed that the QPO properties above 100 keV are predom-inantly influenced by the hard tail component, likely originating from a relativistic jet.Since the high-energy QPO waveform can be well constructed using the phase function derived from the low-energy light curves, a jet-corona coupling precession scenario could explain our findings.Fitting the phase-folded profile with jet precession model suggested a jet velocity of ∼ 0.4c.
We are grateful to the anonymous referee for constructive comments that helped us improve this paper.This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA's Goddard Space Flight Center, and the Insight-HXMT mission, a project funded by China National Space Administration (CNSA) and the Chinese Academy of Sciences (CAS).This work is supported by the National Key R&D Program of China (2021YFA0718500) and the National Natural Science Foundation of China under grants, U1838201, U1838202, 12173103, U2038101 and U1938103.This work is partially supported by International Partnership Program of Chinese Academy of Sciences (Grant No.113111KYSB20190020). Ling D. Kong is grateful for the financial support provided bu the Sino-German (CSC-DAAD) Postdoc Scholarship Program (57251553).VMD is a recently proposed technique designed to decompose a time series into several intrinsic mode functions (Dragomiretskiy & Zosso 2014).In the HHT method, the original mode decomposition technique is the EMD (Huang et al. 1998), which is an algorithm based on interpolation of local extrema.However, this algorithm has a significant drawback: the problem of mode mixing, where widely different scales could appear in a single IMF component, or a coherent signal fragmented into separate parts that appear in more than one IMF component.In comparison with the traditional EMD method, the VMD algorithm theoretically reduces the impact of mode mixing by simultaneously decomposing the input signal into a sum of IMFs while analytically determining the limited center frequency and bandwidth.
In the VMD, the IMFs are represented as amplitude-modulated-frequency-modulated (AM-FM) signals, given by u k (t) = A k (t) cos (ϕ k (t)), where the phase for the kth mode, ϕ k (t), is a non-decreasing function, and the envelope A k (t) is non-negative (A k (t) ≥ 0).When applied to an input signal f (t), the VMD initializes the set of all modes as {u k } = {u 1 , ..., u K } along with their corresponding center frequencies {ω k } = {ω 1 , ..., ω K }.The goal is to minimize the sum of bandwidths of IMFs by solving the optimization problem (1).To estimate the bandwidth of each mode (BK k ), the VMD employs several signal processing tools.Details on these signal processing tools can be found in Dragomiretskiy & Zosso (2014).The bandwidth of each mode is estimated using H 1 Gaussian smoothness of the shifted signal, i.e. the L 2 norm of the gradient: where δ is the Dirac function, and ' * ' denotes convolution.The term enclosed in the square brackets represents the analytic signal of u k (t).With Lagrangian multipliers, λ, the augmented Lagrangian, L, can be expressed as The solution to the original minimization problem ( 1) is now achieved as the saddle point of the augmented Lagrangian, which can be found using the alternate direction multiplier method (ADMM).Then all the modes gained from solutions in Fourier domain are written as where f (ω) and ûk (ω) represent the Fourier transform of the signal f (t) and mode u k (t), respectively, ω is the Fourier frequency.The VMD algorithm iteratively refines the modes {u k } and their corresponding center frequencies {ω k } until they stabilize.In this algorithm, the center frequency ω k is determined within the Fourier domain by the formula: This approach effectively extends Wiener filtering to accommodate adaptive filtering across multiple frequency bands.For a comprehensive exploration of this technique and the constrained variational optimization problem, see Dragomiretskiy & Zosso (2014).

B. SIMULATION ANALYSIS
To assess the reliability of the QPO phase-folding analysis results, we conduct a comprehensive simulation analysis for our HHT-based method to test its robustness.This simulation analysis consists of three parts: (1) testing the performance of the VMD algorithm to verify that a simulated QPO, as strong as the observational one, is indeed decomposed accurately; (2) exploring how the VMD algorithm deals with a "pure noise" signal; and (3) addressing a potential issue that may arise from (1) and (2) in the phase-folding analysis.
While it has now become standard practice to fit power spectra in X-ray binaries with a combination of both broad and peaked components modeled as Lorentzians (see e.g.Belloni et al. 2002;Yang et al. 2022;Zhou et al. 2022;Shui et al. 2023b), there is currently no physical backing behind the use of Lorentzian components.Additionally, we find the continuum of the power spectrum from obsID P011453500201 in 30-120 keV can be well fitted with a broken power law function (see Figure 7).Consequently, we individually model the red noise with a broken power law function and QPO harmonics with the Lorentzian functions.Using the function shapes derived from the modeling of the observed power spectrum, we apply the algorithm from Timmer & Koenig (1995) to generate 1684-s light curves with a 0.05-s time resolution for both red noise and QPO components.In this algorithm, the power spectral amplitude at each frequency is randomly drawn from a χ 2 distribution, and not fixed at the amplitude of the underlying power spectrum.Additionally, the Fourier phases given to the signal by this algorithm is random.For the broken power law funtion of the red noise spectrum, we set the two power law indices at 0.28 and 1.59, respectively, the break frequency at 5.61 Hz, and the fractional rms amplitude at 30%.As for the Lorentzian function of the QPO component, we set the centroid frequency at 1.78 Hz, the FWHM at 0.12 Hz, and the fractional rms amplitude at 15%, respectively.The simulations of the light curves are performed using Stingray5 , which is an open-source Python package (Huppenkothen et al. 2019).In this study, we focus on the QPO fundamental and therefore do not include the harmonic component in the simulation analysis.The total light curve with an average count rate of 600 cts s −1 consists of the red noise, QPO and Poisson noise, which is presented in the left panel of Figure 8. Subsequently, the VMD technique is used to decompose this light curve.It is also important to note that Poisson noise is the observational errors caused by statistical fluctuations in photon counting, while red noise is related to the intrinsic stochastic behavior of the source.Therefore, Poisson noise is specifically incorporated by Poisson sampling of the variable source flux which includes both the red noise and QPO components.The right panel of Figure 8 shows the power spectra of the simulated total (in black), QPO (in grey) light curves and decomposed QPO IMF (in blue).As one can see, the power spectrum of the decomposed QPO IMF is well consistent with that of the simulated QPO signal.This demonstrates that the VMD can precisely recover the essential features of the QPO in the Fourier domain, including the centroid frequency, width and amplitude.We then perform additional Hilbert transform analyses on both the decomposed and simulated QPO light curves.This is followed by comparing the trends in QPO count rate and instantaneous phase between the two sets of light curves.The top panels of Figure 9  and bottom panels represent the case that the two signals are in perfect alignment.We find the structures of the decomposed QPO are highly consistent with those of the simulated QPO signal, but with local fluctuations, possibly due to mode mixing between the QPO and noise components, or the statistical fluctuations.Specifically, in the right middle panel of Figure 9, a small number of data points are distributed in the upper left and lower right of the plot, corresponding to the phase difference distributed around −2π and 2π in the right bottom panel.Since the phase of the QPO varies from −π to π, this phenomenon reflects that the phases of the simulated and decomposed QPO signals in a small part of time bins have local shifts.
The VMD technique is also applied to a pure Poisson noise signal.In Figure 10, we plot the power spectra of the Poisson noise and corresponding VMD IMFs.The result shows that the decomposed IMFs are evenly distributed in frequencies, suggesting that VMD is essentially a set of adaptive Wiener filter bank capable of separating modes with different center frequencies.The equivalent filter bank properties of the VMD have been extensively studied in literature (see e.g.Wang et al. 2015).From another perspective, the VMD algorithm is capable of generating the output IMFs for any given signal, even if the signal is pure noise.This indicates that the IMF corresponding to the QPO component might be influenced by the noises, potentially leading to local variations between the real and decomposed QPO signals observed in Figure 9.To explore the potential impact of the local variations between simulated and decomposed QPO light curves presented in Figure 9 on the phase-folded light curve, we conduct the phase-folding on the initially simulated total and QPOsubtracted light curves, utilizing the instantaneous phase of the decomposed QPO IMF.Here, the QPO-subtracted light curve refers to the red noise with its Poisson noise.In our observational study, we combined data from seven different observations, each with varying QPO frequencies and exposures.Our simulation analysis is designed to closely replicate the observational analysis, taking into account these variations in QPO frequencies and exposures as accurately as possible.Therefore, we simulate seven sets of light curves in which the QPO frequency and exposure times are set to be the same as those in the corresponding observations.For each total light curve, we apply the VMD method and the Hilbert transform.Subsequently, the phase-folded total and QPO-subtracted light curves are obtained by phase-folding across the seven light curves using the QPO instantaneous phases obtained from the HHT analysis.As shown in Figure 11, the phase-folding profile of the total light curve exhibits a significant modulation, which is defined as the QPO waveform.The amplitude of the modulation is estimated to be ∼ 0.21, with a corresponding rms of ∼ 14.9%, consistent with the initially set value in the simulation (15%).However, a relatively weaker modulation (with an amplitude of ∼ 0.04) of phase-folded QPO-subtracted light curve is also observed.This indicates that a small portion of the red noise is decomposed into the QPO IMF, resulting from the issue of mode mixing.
As the phase-folded noise light curve simulated in the low-energy band can exhibit modulations due to mode mixing (see Figure 11), it is important to confirm that the observed modulations in the high-energy bands (see Figure 3) are indeed from the QPO component and not from the red noise.To achieve this, we simulate high-energy light curves for the seven observations based on the conditions of the observational data.Specifically, in the 150-170 and 170-200 keV energy bands, the average total count rates of the observed light curves are ∼ 32 and 47 cts s −1 , respectively, with the net count rates of ∼ 2 cts s −1 in both energy bands.For simplicity, we assume that the light curves in different energy bands for both the QPO and red noise are fully coherent.This indicates that the light-curves of both the QPO and red noise have consistent phase functions across different energy bands.Additionally, we assume that the fractional rms of the red noise and the QPO remains constant with photon energy.Therefore, to simulate red noise and QPO light curves in the two high-energy bands, we scale the absolute amplitude of the original simulated red noise and QPO light curves in the 30-120 keV energy band by a factor of 1/300 (∼ 2 cts s −1 /∼ 600 cts s −1 ).Additionally, in these simulated total light curves, Poisson noise is also incorporated.
Following the analysis presented in Section 3, for each energy band, we combine these seven light curves to generate the PDS.Similar to the observational findings presented in Figure 4, the combined power spectra of the simulated light curves in the high-energy bands are dominant with Poisson noise, resulting in QPO peaks that are too weak to be robustly detected (see Figure 12).To rigorously evaluate the phase-folding profiles under the impact of Poisson noise, we independently simulate 10 4 light curves for each of the two energy bands.Then, the phase-folding is performed on the simulated total and QPO-subtracted light curves for both energy bands.It should be noted that, due to the low SNR in the high-energy band, conducting the HHT analysis there is not feasible.Consequently, the phase-folding of the high-energy light curves is carried out using the instantaneous phase obtained from the HHT analysis in the low-energy band.
The upper panels of Figure 13 present phase-folding profiles using 1000 samples from the simulated light curves in the two energy bands.It is evident that in both energy bands, the phase-folding profiles of the total light curves exhibit significant modulations, while those of the QPO-subtracted light curves are relatively weaker, likely being dominated by the random fluctuations of Poisson noise.Based on the cross-correlation technique described in Appendix C, for each sample, we calculate the maximum cross-correlations between phase-folding profiles of the total, QPO-subtracted, pure Poisson noise light curves in the high-energy bands and the reference profile.The reference profile is the phase-folded light curve in the low-energy band, plotted in blue in Figure 11.
As shown in the bottom panels of Figure 13, the distributions of maximum cross-correlations of the pure Poisson noise and QPO-subtracted signals exhibit significant overlap in both energy bands.This demonstrates the QPO-subtracted signals in the two energy bands are dominated by the Poisson noise.Therefore, the modulations in the phase-folding profiles of the QPO-subtracted signals in the high-energy bands could only be determined at the confidence level of < 2σ, which is inconsistent with our observational results (∼6-8 σ).Moreover, in a real case, the red noise light curves from two energy bands are not perfectly correlated, indicating that the modulations of the QPO-subtracted signals in the high energy bands could be even weaker than the simulated ones.However, the distributions of the total signals (with the QPO component) are significantly different from those of the pure Poisson and QPO-subtracted signals, indicating that the modulations of total light curves in the high energy bands could be detected at the ∼5-8 σ confidence level, consistent with the observational results.
In the simulation analysis above, we assume that the fractional rms of QPO and red noise does not vary with photon energy.However, in reality, the fractional rms of both signals may indeed change with photon energy.Therefore, we conduct tests to examine how the amplitude of the QPO and red noise affects the phase-folding method performed in the 150-170 keV and 170-200 keV energy bands.Firstly, we keep the fractional rms of the red noise constant at the value observed in the 30-120 keV energy band, while varying the QPO fractional rms.The top two panels of Figure 14 present the dependence of the modulation significance on the QPO rms.The significance is determined with the Monte Carlo method described in Appendix C. It is evident that the modulation significance of the total profiles increases with the QPO rms, whereas that of the QPO-subtracted profiles remain relatively constant (0-4σ).When the QPO fractional rms is lower than ∼ 8%, the modulation of the total profiles cannot be robustly determined (< 5σ).Subsequently, we fix the fractional rms of the QPO signal at the value of that in the 30-120 keV energy band, and vary the fractional rms of the red noise.As shown in the two bottom panels of Figure 14, both the total and QPO-subtracted profiles exhibit modulations that can be determined with higher confidence levels when the fractional rms of the red noise is higher.However, for the cases that the fractional rms of the red noise is lower than 60%, the modulations of the QPO-subtracted profiles could only be determined at the confidence level of < 5σ, which is inconsistent with our observational results (∼6-8 σ).Moreover, previous studies on the red noise in BHXRBs have indicated that its fractional rms generally decreases with photon energy (see e.g.Bu et al. 2021;Yang et al. 2022).This suggests that the fractional rms of the red noise is likely to be lower than 30% in the 150-170 keV and 170-200 keV energy bands, and could not significantly influence the determination of the intrinsic QPO modulation.
In our simulation analysis, we model the red noise component using a broken power law function.We have also conducted the simulation analysis using the zero-frequency Lorentzian function to model the red noise, and obtained similar results.Our simulation analyses provide supports for the idea that the observed modulations of the phasefolding profile in the energy ranges of 150-170 and 170-200 keV could indeed originate from the QPO component, despite the presence of some mode mixing issues in the VMD.

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

1Figure 1 .
Figure 1.HE light curves of MAXI J1535-571.In this plot, red triangles denote observations with type-C QPO detection, while blue circles represent those without type-C QPO detection.Each point represents one observation, and the shaded area indicates the selected observations for detailed timing analyses.

Figure 2 .
Figure2.Hilbert-Huang transform analysis of a representative example of a 14-s-long light curve with a 0.05s time resolution of HE in energy range of 30-120 keV.a, representative example of a 14-s-long light curve and corresponding five IMFs obtained from VMD algorithm.b, power density spectra of the original HE light curve, as well as the second and third IMFs.c, the QPO (the second) IMF and corresponding phase function, instantaneous frequency and amplitude.

Figure 3 .
Figure 3.The Hilbert-Huang transform (HHT) phasefolded QPO waveforms of MAXI J1535-571 from Insight-HXMT/HE light curves.To enhance the SNR, waveforms from all selected observations have been collectively combined, and the final combined waveform are normalized by the mean value of the count rate.These QPO phases are obtained through HHT analysis of light curves in the 30-120 keV energy band.Each panel showcases the significance level determined via the cross-correlation technique, as detailed in Appendix C.

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

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

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

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

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

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

Figure 12 .
Figure 12.Power density spectra of the simulated light curves in the energy ranges of 150-170 (in orange) and 170-200 (in blue) keV.

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

Table 1 .
Log of Insight-HXMT Observations Used in This Work Russell et al. (2019)0,150-170and 170-200keV energy bands, respectively.The phase differences between these energy bands represent soft lags.Based on the radio observations,Russell et al. (2019)suggested an inclination of ⩽ 45 • .We also perform fittings with the jet precession model, assuming i = 45 • .