Detection of a Low-frequency Quasi-periodic Oscillation in the Soft State of Cygnus X-1 with Insight-HXMT

Cygnus X-1 is a bright persistent X-ray binary consisting of a black hole accreting matter from an OB star (Bowyer et al. 1965; Bolton 1972; Webster & Murdin 1972). Since its discovery, Cygnus X-1 has been the target of extensive multiwavelength observational campaigns (e.g., Herrero et al. 1995; Nowak et al. 1999; Gallo et al. 2005; Albert et al. 2007; Grinberg et al. 2013; Kantzas et al. 2021), making it one of the most-studied X-ray sources. From the very early observations of Cygnus X-1, it was clear that the X-ray emission undergoes dramatic changes in spectral distribution and brightness (Tananbaum et al. 1972).

Analogous changes have been observed in most later discovered black hole X-ray binaries (BH XRBs) and have been classified into different accretion states (e.g., Remillard & McClintock 2006; Belloni 2010). The accretion states that have been first identified are the low-hard state (LHS) and the high-soft state (HSS). In the LHS, the energy spectrum can be described by a power law (the photon index ∼1.5–2.1) and in the HSS by a multitemperature blackbody component. It is generally accepted that the multitemperature blackbody emission arises from an accretion disk around the black hole (Shakura & Sunyaev 1973) and the power law is the result of Compton upscattering of cool photons by hot electrons close to the black hole (see reviews in Done et al. 2007; Yuan & Narayan 2014). Accretion states are also identified according to the timing properties of the source (see reviews by Remillard & McClintock 2006; Belloni 2010; Belloni & Stella 2014). The power density spectrum (PDS) in the LHS usually exhibits low-frequency quasi-periodic oscillations (LF QPOs, ∼0.01–10 Hz) on top of broadband noise. QPOs are classified as types A, B, and C according to their frequency, amplitude, and coherence (see Wijnands et al. 1999; Casella et al. 2004 for detailed definitions of QPO types). Strong high-coherence type-C QPOs are usually observed in the LHS, when the total fractional rms amplitude is high (∼30%). In the HSS the amplitude of the variability drops dramatically (≲1%) and no QPOs are observed. The above phenomenon has been observed in almost every transient BH XRB, which has the largest population of known Galactic BHs. The LHS is usually observed at the beginning and the end of periods of enhanced accretion (outbursts), when the luminosity is low, while the HSS occurs at maximum luminosity. Between these two extremes, several transition states are also observed (Remillard & McClintock 2006; Belloni 2010).

However, some persistent BH XRBs can exhibit accretion states different from the transient BH XRBs. Cygnus X-1 is one of these sources, as its HSS is characterized by both multitemperature blackbody and Compton power-law emission (e.g., Gierliński et al. 1999). The X-ray timing properties of accretion states in Cygnus X-1 differ from the general description mentioned above. Type-C QPOs have never been detected during the LHS of Cygnus X-1 (Grinberg et al. 2014; Ingram & Motta 2019) and the HSS shows high variability (fractional rms amplitude ∼25%; Axelsson et al. 2005; Grinberg et al. 2014). However, there were broad and weak narrow humps detected in its PDS during different spectral states, which were sometimes called QPOs (e.g., Cui et al. 1997b; Paul et al. 1998; Pottschmidt et al. 2003; Axelsson et al. 2005; Shaposhnikov & Titarchuk 2006; Grinberg et al. 2014).

In this paper, we report the detection of a millihertz QPO in one of the Insight-HXMT exposures of Cygnus X-1 during the soft state. In Section 2, we describe the data reduction, in Section 3 the spectral and timing analysis, and in Section 4 we compare the detected signal to similar signals observed in other BH XRBs and discuss possible physical scenarios.

The Insight-Hard X-ray Modulation Telescope (HXMT) (Zhang et al. 2020) is China's first X-ray astronomy satellite. Insight-HXMT has three major scientific payloads: the Low Energy Telescope (LE, 1–10 keV), the Medium Energy Telescope (ME, 5–35 keV), and the High Energy Telescope (HE, 20–250 keV), capable of observations down to a time resolution of 1 ms, 276 μs, and 25 μs, respectively. The effective area of the three payloads is 384 cm², 952 cm², and 5000 cm², respectively.

In order to reduce the data file size, a long Insight-HXMT observation is man made, split into multiple segments. Each segment lasts approximately three hours, which is named "exposure" and identified with an Exposure ID. ¹ We analyzed the longest observation of Cygnus X-1 (Obs. ID P0101315001, start time 2017 August 24T02:49:21 UTC) in the Insight-HXMT data archive. The selected observation contains 13 exposures, for a total duration of about 150 ks (see Figure 1).

Zoom In Zoom Out Reset image size

We performed data reduction by using the Insight-HXMT Data Analysis Software package (HXMTDAS) V2.02. For all three payloads, we extracted cleaned events applying the following filtering criteria: (1) exclusion of time periods when Insight-HXMT passes through the South Atlantic Anomaly (T_SAA > 300 && TN_SAA > 300 && SAA_FLAG==0); (2) elevation angle >10° (ELV > 10); (3) pointing offset angle <0.1° (ANG_DIST < 0.1); and (4) cutoff rigidity > 8 GeV (COR > 8). An additional filtering criterion was applied for LE data: elevation angle from Earth bright limb >20° (DYE_ELV > 20°). For LE, ME, and HE, we computed binned light curves, energy spectra, and corresponding background from cleaned event files using the HXMTDAS tools <le/me/he > lcgen, <le/me/he > specgen, and <le/me/he > bkgmap, respectively. In particular, we extracted source and background light curves with time resolution 1/256 (≈0.0039) s in the energy bands 1–10 keV, 5–30 keV, and 20–250 keV for LE, ME, and HE, respectively.

We also downloaded the daily light curves from the X-ray all-sky monitoring instruments Swift/BAT (Krimm et al. 2013) and MAXI (Matsuoka et al. 2009), selecting Swift/BAT and MAXI data simultaneous to our Insight-HXMT exposure (MJD 57989L; see Figure 2).

3.1. Timing Analysis

For each payload and exposure, we subtracted the background from the light curve, we divided the light curves into segments of 128 s, and, for each segment, we computed Leahy-normalized PDSs (Leahy et al. 1983) using the HEASARC tool powerspec (XRONOS package). PDSs were averaged, obtaining in this way a single PDS for each exposure. Because of the chosen time resolution and segment length, PDSs have frequency resolution d ν = 1/128 = 0.0078125 Hz and Nyquist frequency ν_N = 128 Hz.

Zoom In Zoom Out Reset image size

Figure 1 shows the light curve with a time resolution of 100 s for each payload for the entire selected observation. We visually inspected PDSs for each exposure and find obvious QPO features in all three payloads during the first exposure (Exp. ID P010131500101-20170824-01-01; red region in Figure 1). In this exposure, the PDS of each payload is characterized by a power-law broadband noise with a narrow feature on top of it at low frequency (Figure 3). This narrow feature appears to be more prominent at higher energies (ME and HE), and it is not present in any other exposures of the selected observation. To ascertain the presence of the narrow feature in the LE low-energy regime (1–10 kev), we performed the same data reduction and analysis described above in two sub-low-energy bands, 1–3 keV and 3–10 keV, respectively.

Zoom In Zoom Out Reset image size

We then subtracted the contribution from the Poisson noise, which was estimated by averaging Leahy-normalized PDSs above 100 Hz. The noise levels (1.99, 2.00, 1.99, and 2.09 for LE 1–3 keV, LE 3–10 keV, ME, and HE) are consistent with the expected Poisson noise level in Leahy-normalized PDSs (Leahy et al. 1983). The average net count rates of 870.50,178.00, 87.04, and 52.60 counts/s for LE(1–3 keV), LE(3–10 keV), ME, and HE are used for calculating rms-normalized PDSs (Belloni & Hasinger 1990). The total fractional rms amplitude of each PSD is 11.4%, 27.3%, 28.7%, and 37.1% (below 32 Hz). We then fit the PDSs using a model consisting of a power law plus a narrow Lorentzian. The fitting was performed in XSPEC (Arnaud 1996). We obtained a statistically acceptable fit for the HE PDS (see Figure 3 and Table 1), while we noticed broad frequency residuals around 0.5 Hz for the LE and ME PDS. Adding a broad Lorentzian component (BLN; a Lorentzian function with ν₀ fixed at 0) in the model, we obtain a statistically acceptable fit also for LE and ME (see Figure 3 and Table 1). We then calculate the ${\nu }_{\max }$ ($=\sqrt{{\nu }_{0}^{2}+{(\mathrm{FWHM}/2)}^{2}}$) to represent its characteristic frequency (e.g., Belloni et al. 2002), the values of which are ${0.50}_{-0.07}^{+0.07}$, ${0.56}_{-0.07}^{+0.08}$, and ${0.56}_{-0.11}^{+0.13}$ for LE (1–3 keV), LE (3–10 keV), and ME data.

Table 1. Best-fitting Parameters of the LE and ME PDSs

Component	Parameter	LE	LE	ME	HE
		(1–3 keV)	(3–10 keV)
POW	Index	${1.16}_{-0.04}^{+0.04}$	${1.08}_{-0.03}^{+0.03}$	${1.08}_{-0.04}^{+0.04}$	${0.95}_{-0.02}^{+0.02}$
LOR	ν0 (mHz)	${89.27}_{-2.87}^{+5.68}$	${86.91}_{-1.71}^{+1.47}$	${87.53}_{-2.63}^{+1.70}$	${88.84}_{-1.48}^{+0.72}$
(QPO)	FWHM (mHz)	${14.03}_{-9.35}^{+15.53}$	${7.39}_{-3.75}^{+4.86}$	${2.88}_{-2.88}^{+3.62}$	${2.64}_{-2.64}^{+4.19}$
	rms(%)	${2.18}_{-0.47}^{+0.44}$	${5.82}_{-0.87}^{+0.76}$	${4.75}_{-0.91}^{+0.79}$	${5.94}_{-0.94}^{+0.89}$w
LOR	ν0 (Hz)	0	0	0	⋯
(BLN)	FWHM (Hz)	${0.99}_{-0.13}^{+0.15}$	${1.12}_{-0.13}^{+0.15}$	${1.11}_{-0.22}^{+0.26}$	⋯
	rms(%)	${5.57}_{-0.41}^{+0.39}$	${14.31}_{-0.77}^{+0.74}$	${13.39}_{-1.18}^{+1.12}$	⋯
	χ2/dof	161.08/156	189.59/156	163.90/156	178.93/158

Note. The model is Powerlaw+Lorentzian+Lorentzian for LE and ME, Powerlaw+Lorentzian for HE. The uncertainties are the 1σ ranges.

Download table as:  ASCIITypeset image

All the best-fitting parameters are listed in Table 1. The power-law noise with an index of ∼1 and the fractional rms are consistent with the soft state of Cygnus X-1 (e.g., Cui et al. 1997a; Axelsson et al. 2005, 2006; Grinberg et al. 2014). A narrow QPO is detected in all PDSs with the significance of 4.7σ, 6.8σ, 5.6σ, and 6.3σ for the LE (1–3 keV), LE (3–10 keV), ME, and HE, respectively. The QPO frequencies at different energy bands are consistent within uncertainties; the average value is ∼88 mHz. The QPO is very narrow, especially at higher-energy bands (see the FWHM values in Table 1). The quality factors (Q = ν₀/FWHM) are 10.8 ± 7.7, 15.3 ± 8.1, >13.7, and >13 for the LE(1–3 keV), LE(3–10 keV), ME, and HE, respectively. We further investigated the energy dependence of the QPO rms with more subenergy bands. We first produced the PDSs in the energy bands of 1–3 keV, 3–5 keV, and 5–10 keV for LE data; 5–10 keV, 10–15 keV, and 15–30 keV for ME data; and 30–50 keV and 50–150 keV for HE data. We then calculated the fractional rms of the QPO by fitting each PDS. The QPO rms becomes almost constant (∼5%) above 3 keV and up to at least 50 keV (see Figure 4).

Zoom In Zoom Out Reset image size

3.2. Further Assessing the Significance of the QPO

The classical method to determine the detection level of a signal is based on the statistical properties of the white noise (van der Klis 1989). For an averaged PDS, the white-noise power follows a χ² distribution with 2MW dof, where M is the number of averaged powers at a specific frequency and W is the number of averaged frequency bins. Integrating the χ² probability density function, it is possible to determine the p-value of a given power and define detection levels above the underlying noise. However, this method is not reliable for accessing detection above noise processes differing from white noise (Israel & Stella 1996; Vaughan 2005, 2010).

In our case, the QPO is detected in a frequency range dominated by red noise. To access the significance of the QPO with respect to the red noise, we applied the method proposed by Vaughan (2010) and implemented in the Python package stingray (Huppenkothen et al. 2019). We only considered PDS below 8 Hz, where the red noise dominates. We fitted the PDSs with a power-law function plus a constant and sample the model parameters with a Markov Chain Monte Carlo method in the emcee package (Foreman-Mackey et al. 2013). We then generated 10⁵ fake PDSs using the sampled parameters, in order to derive a posterior probability distribution for the test statistics and thus a posterior predictive p-value. Our observed PDSs are averaged by M 128s segments, where M is 18, 22, and 27 for LE, ME, and HE. We considered the M when generating the fake PDSs and did not rebin the fake PDSs. So, the derived p-value has been corrected for the number of periodograms searched including the number of segments and the number of frequencies.

We fitted each simulated PDS with a power-law plus a constant model and calculate (2I_j /S_j ), where I_j is the power at a given frequency f_j of a PDS and the S_j is the power from the model. We then plot the 99.450%, 99.730%, and 99.993% (2σ, 3σ, and 4σ) significance levels, which is obtained from 10⁵ samples of (2I_j /S_j ) at each frequency (Figure 5). There is only one observed power at the frequency ∼86 mHz that is larger than that of the 99.993% simulated PDS, which indicates a significant signal (>4σ) around this frequency in the PDS of LE (3–10keV), ME, and HE. In the PDS of LE (1–3 keV), the power at ∼86 mHz is marginally at 2σ.

Zoom In Zoom Out Reset image size

We first defined the test statistic ${T}_{{\rm{R}}}=\max (2{I}_{j}/{S}_{j})$. Then the p-values of the observed ${T}_{{\rm{R}}}^{\mathrm{obs}}$ from the posterior probability distribution of T_R are 0.00394, 0.00334, and 0.00027 for LE (3–10 keV), ME, and HE, respectively, in Figure 6. The combined p-value of the three is 7.45 × 10⁻⁷ by using the Fisher method as they are independent tests for the three instruments. The frequency with maximum T_R from the three instruments are the same, so the p-value at this specific frequency should be smaller than the combined p-value (e.g., Huppenkothen et al. 2017), demonstrating a very significant signal of ∼86 mHz existing above the red noise.

Zoom In Zoom Out Reset image size

We also used the likelihood ratio test (LRT) to compare the two models with/and without a Lorentzian component. The likelihood ratio is defined as ${T}_{\mathrm{LRT}}=-2\mathrm{log}({L}_{0}/{L}_{1})$, where the L₀ is the likelihood of the simple model and the L₁ is the likelihood of the complex model. We fitted each simulated PDS with the two models and computed the T_LRT. The posterior probability distributions of the T_LRT can give us the p-values of ${T}_{\mathrm{LRT}}^{\mathrm{obs}}$ of 0.00008, 0.00016, and 0.00001 for LE (3–10 keV), ME, and HE, respectively (Figure 6). Therefore, the statistics of T_R and T_LRT both demonstrate a significant QPO detection in LE (3–10 keV), ME, and HE. The p-values of T_R and T_LRT of LE (1–3 keV) are at least one order of magnitude larger than those of higher-energy bands, which means the QPO is less significant below 3 keV.

3.3. Spectral and Color Analysis

We performed spectral analysis of the selected exposure by combining the data from LE, ME, and HE. The effective areas of the three instruments have been calibrated using Crab nebula observations as described in Li et al. (2020). We jointly fit energy spectra of LE (1–10 keV), ME (8–30 keV), and HE (30–150 keV) with the model const∗tbabs∗diskbb+reflkerr) (see Figure 7), where abundances and cross sections of the absorption by the Galactic interstellar medium are set according to Wilms et al. (2000) and Verner et al. (1996). reflkerr is a relativistic reflection model that computes direct and reflection spectra assuming either a slab or a spherical plasma geometry (Niedźwiecki et al. 2019). In this framework, the Comptonization spectrum is computed with compps (Poutanen & Svensson 1996). We fix the black hole spin (0.99, Zhao et al. 2020, 2021), the inclination angle (275), Miller-Jones et al. 2021), and the index of the outer disk q_out (3), and we leave the inner disk index q_in and break radius R_br free (Fabian et al. 2012; Walton et al. 2016). We also tied the temperature of seed photons for Comptonization with the disk temperature. All of the best-fitting parameters are listed in Table 2. According to the electron temperature kT_e and optical depth τ of the corona, we can derive the photon index of the Comptonization component (Zdziarski et al. 1996). The derived photon index Γ = 2.74 ± 0.11 and the best-fitting disk temperature T_in = 0.49 ± 0.01 keV are consistent with the soft state of Cygnus X-1 (Tomsick et al. 2014; Walton et al. 2016; Kawano et al. 2017; Lubiński et al. 2020). The unabsorbed X-ray flux (0.1–100 keV) is 8.16 × 10⁻⁸ erg s⁻¹ cm⁻², which corresponds to ∼2% Eddington luminosity at a distance of 2.22 kpc and a BH mass of 21.2 M_⊙ (Miller-Jones et al. 2021).

Zoom In Zoom Out Reset image size

Table 2. Best-fitting Parameters of the Energy Spectrum

Component	Parameter		Value
TBABS	nH	[1022 cm−1]	${0.46}_{-0.04}^{+0.04}$
DISKBB	kTin	[keV]	0.49${}_{-0.01}^{+0.01}$
	Norm	[104]	${2.64}_{-0.19}^{+0.17}$
REFLKERR	qin		${7.52}_{-1.75}^{+0.58}$
	Rbr	[rg ]	${3.29}_{-0.21}^{+0.23}$
	Rin	[rg ]	${2.22}_{-0.16}^{+0.10}$
	τ		${0.59}_{-0.03}^{+0.04}$
	Afe	[solar]	<10
	kTe	[keV]	${87.74}_{-2.39}^{+2.82}$
	$\mathrm{log}\xi$	log[erg cm s−1]	${4.33}_{-0.03}^{+0.06}$
	frefl		${0.44}_{-0.07}^{+0.08}$
χ2/dof			1332.78/144

Note. The uncertainties are the 1σ ranges.

Download table as:  ASCIITypeset image

To confirm that Cygnus X-1 is in the soft state in our selected exposure, we computed count rates from Swift/BAT and MAXI light curves simultaneous to our data (MJD 57989). The count rate in the 15–50 keV (BAT) and 2–4 keV (MAXI) energy band is (4.73 ± 0.28) × 10⁻² counts cm⁻² s⁻¹ and 1.97 ± 0.23 counts s⁻¹, respectively. The corresponding hardness ratio (4–10 keV/2–4 keV) is 0.39 ± 0.08. According to the state classification criteria of Cygnus X-1 described in Grinberg et al. (2013), the observed X-ray intensity and hardness ratio are both consistent with the source being in the soft state.

We detected a very low-frequency (∼88 mHz) QPO only during the first 12 ks of the longest Insight-HXMT observation of Cygnus X-1. This QPO is significantly detected in LE, ME, and HE independently. The fractional rms amplitude is energy independent (above 3 keV; see Figure 4). The broadband noise at a low frequency of the PDSs from all three Insight-HXMT instruments is dominated by a power-law component with index ∼1 (Table 1), which is consistent with the soft state of Cygnus X-1 (e.g., Axelsson et al. 2005). The measured fractional rms amplitude above 10 keV is more than two times larger than that below 3 keV (see Section 3.1), which is also a characteristic of the soft state of Cygnus X-1 (e.g., Grinberg et al. 2014). Both the X-ray spectral properties of and all-sky monitoring data corresponding to this exposure also show that Cygnus X-1 is in the soft state (Table 2 and Figure 2). About state classification, it is worth mentioning that Cygnus X-1 has never reached the so-called thermal-dominated or soft state as defined according to BH low-mass X-ray binary phenomenology (BH LMXBs; e.g., Remillard & McClintock 2006; Belloni 2010). The X-ray spectra of the soft states of other BH LMXBs are dominated by a thermal disk component, while in Cygnus X-1 strong Comptonization and reflection components are also significant (Walton et al. 2016; Kawano et al. 2017; Lubiński et al. 2020). Furthermore, the fractional rms amplitude measured from soft-state PDSs is usually ∼20%–30% in Cygnus X-1 (e.g., Axelsson et al. 2006; Grinberg et al. 2014), much larger than the one measured in BH LMXBs in the thermal-dominated state. Keeping in mind these differences, we will refer to this state of Cygnus X-1 simply as the soft state.

Low-frequency QPOs are quite common in the hard state of BH XRBs. In particular, type-C QPOs are observed with fractional rms amplitude up to 20% and quality factor Q (ν₀/FWHM) ≥ 10 (Casella et al. 2005; Motta 2016). Type-C QPOs are sometimes also detected in the soft state (e.g., XTE J1550−564, GRO J1655−40), but with a higher characteristic frequency (≥20 Hz) and smaller rms amplitude than those seen in the hard state (Homan et al. 2001; Motta et al. 2012). Type-C QPOs with similar characteristics are systematically observed in most of the known BH LMXBs. However, type-C QPOs have never been detected in Cygnus X-1 (Belloni 2010; Ingram & Motta 2019). Pottschmidt et al. (2003) detected weak narrow QPOs (ν₀ ∼ 0.1–1 Hz) sometimes during the soft state. Rapisarda et al. (2017) also detected a millihertz QPO (ν₀ ∼ 60 mHz) during the soft state. The millihertz QPO we detected has a large quality factor Q (>13; see ν₀ and FWHM in Table 1), which is similar to that of type-C QPOs. However, its frequency and rms are smaller than those of most type-C QPOs observed in other BH XRBs (e.g., Motta et al. 2015). On the other hand, the rms of type-C QPOs usually increases with photon energy and becomes constant above tens of keV (e.g., Casella et al. 2004; Huang et al. 2018), while the rms of the millihertz QPO we detected in Cygnus X-1 does not show an energy dependence above 3 keV (see Figure 4 and Table 1). Another important difference between type-C QPOs and the millihertz QPO we detected is that the former can be observed for relatively long periods of time (∼days to tens of days), while the latter appears only in the first 12 ks of one observation. These fundamental differences suggest that the physical mechanism responsible for the observed millihertz QPO in Cygnus X-1 is different from the one producing type-C QPOs in BH LMXBs.

The origin of millihertz QPOs (<100 mHz) is still not well understood, although they have been detected in different accreting systems, such as accretion-powered X-ray pulsars (e.g., Psaltis 2006) and even ultraluminous X-ray pulsars (e.g., Feng et al. 2010), ultraluminous X-ray sources (e.g., Atapin et al. 2019), and BH XRBs (see Table 2 in Cheng et al. 2019 and references therein). In some BH LMXBs (e.g., H1743−322), millihertz QPOs and type-C QPOs are sometimes detected simultaneously and their characteristics mostly appear to be uncorrelated with each other (Altamirano & Strohmayer 2012). This supports the conclusion that in BH LMXBs, the mechanism producing the millihertz QPO is not related to that producing type-C QPOs. Some of the detected millihertz QPOs appear in high-inclination systems and are thought to have a similar origin to the 1 Hz QPO in dipping neutron star XRBs (e.g., Altamirano & Strohmayer 2012; Armas Padilla et al. 2014). However, the inclination angle of the Cygnus X-1 accretion disk is estimated to be ∼27° (Orosz et al. 2011), making unlikely the relation between our detected millihertz QPO and those observed in high-inclination systems. We have noticed that the millihertz QPOs are all detected during the hard state in BH LMXBs, while they are detected during the soft state of high-mass X-ray binaries (HMXBs; Cheng et al. 2019), which indicates different origins in HMXBs and LMXBs.

The millihertz QPOs in soft states have been previously observed in two other persistent BH HMXBs, namely LMC X-1 (Ebisawa et al. 1989; Alam et al. 2014) and Cyg X-3 (van der Klis & Jansen 1985; Koljonen et al. 2011; Pahari et al. 2017). The millihertz QPOs in Cyg X-3 are thought to be related to major radio flare events (Koljonen et al. 2011), which are quite unique and have not been observed in other BH XRBs, including Cygnus X-1. LMC X-1 share many common spectral and timing properties during the soft state with Cygnus X-1, even though the X-ray Eddington luminosity of Cygnus X-1 is smaller than that of LMC X-1 (e.g., Ruhlen et al. 2011) and LMC X-1 never showed transitions to the hard state. The millihertz QPO in LMC X-1 is also occasionally detected in soft-state observations, and the frequency also varies (Ebisawa et al. 1989; Alam et al. 2014). The rare occurrence in Cygnus X-1, as well as in LMC X-1, favors potential mechanisms that would only happen by chance due to some local inhomogeneity in a system embedded in a wind-accretion environment.

It has been proposed that the millihertz QPOs observed in LMC X-1 are due to global disk oscillations (Titarchuk & Osherovich 2000), i.e., vertical (normal to the disk) disk oscillations triggered by the gravitational force of the central BH. The frequency of such oscillations can be estimated as

$$\begin{eqnarray}&&{f}_{0}\approx 2{\left(\displaystyle \frac{{R}_{\mathrm{in}}}{6{R}_{{\rm{g}}}}\right)}^{-\tfrac{8}{15}}{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}\right)}^{-\tfrac{8}{15}}{\left(\displaystyle \frac{{P}_{\mathrm{orb}}}{3\mathrm{hr}}\right)}^{-\tfrac{7}{15}}\ {\left(\displaystyle \frac{{R}_{\mathrm{adj}}}{{R}_{\mathrm{in}}}\right)}^{-0.3}\,\mathrm{Hz},\end{eqnarray} \tag{ 1 }$$

where R_in is the inner radius of the accretion disk and R_adj is an adjustment radius in the disk, usually of the order of ∼2–3R_in for Shakura–Sunyaev geometrically thin disks (Shakura & Sunyaev 1973; Titarchuk et al. 1998). Using R_adj = 2.5R_in, P_orb = 5.6 days, R_in = 2.2R_g (see Table 2), and M_BH = 21.2M_☉ (Miller-Jones et al. 2021), we obtain an oscillation frequency of ∼86 mHz, roughly consistent with the frequency of the millihertz QPO we detected. In general, the frequency of the global oscillation depends on the size of the disk and, therefore, on both its inner and outer radius. In Equation (1), the disk is assumed to be a geometrically thin Shakura–Sunyaev disk (Shakura & Sunyaev 1973), and the outer disk radius is assumed to be half of the Roche lobe radius (Titarchuk & Osherovich 2000). In this scenario, the detection of similar millihertz QPOs at different frequencies in the same and other sources (Alam et al. 2014; Rapisarda et al. 2017) could be caused by the variation in disk size.

However, as global disk oscillation is the vertical oscillation of the whole disk, the consequently observed modulation is expected to be from the disk emission, where it is mainly in the soft X-ray band (see Figure 7). On the other hand, global disk oscillations would be seen in other BH XRBs not limited to wind accreting systems. The fractional rms amplitude of the detected QPO persisted at a similar amplitude level at least up to 50 keV (see Figure 4), indicating that the QPO either has a hard X-ray spectrum or that the QPO originates primarily in the Comptonization component, such as certain oscillations occur in the disk corona above the standard accretion disk.

We thank Chichuan Jin and Joern Wilms for the helpful discussions. This research has made use of MAXI data provided by RIKEN, JAXA, and the MAXI team. This work was supported in part by the Natural Science Foundation of China (grants 11773055, U1838203, and U1938114), Z.Y. was also supported by the Youth Innovation Promotion Association of CAS (ID 2020265). R.S. acknowledges the support of the PIFI fellowship of CAS under project No. 2019PM0016 and China's Postdoctoral International Exchange Program. W.Y. would like to acknowledge partial support by the National Program on Key Research and Development Project (grant No. 2016YFA0400804).

Facility: Insight-HXMT. -

Software: astropy (Astropy Collaboration et al. 2013), stingray (Huppenkothen et al. 2019), XSPEC (Arnaud 1996).