TRACKING THE EVOLUTION OF QUASI-PERIODIC OSCILLATION IN RE J1034+396 USING THE HILBERT–HUANG TRANSFORM

To obtain the variation of QPO frequency, we applied a modified version of the empirical mode decomposition (EMD), i.e., the ensemble EMD (EEMD) proposed by Wu & Huang (2009), on the combined PN+MOS light curve. The EMD is an iterative "sifting" process for extracting oscillation modes by subtracting the local means from the original data (Huang et al. 1998). The decomposed components are intrinsic mode functions (IMFs) that satisfy the following two conditions: (1) the number of extrema and the number of zero crossings must be identical or differ by one, and (2) the local mean of the data is zero (Huang et al. 1998; Hu et al. 2011). The instantaneous frequency obtained by the Hilbert transform is meaningful only if it is calculated from an IMF. EEMD overcomes the mode-mixing problem, i.e., a modulation with the same timescale distributed across different IMFs, by taking the ensemble mean of the IMFs from combinations of the original data and added white noises (Wu & Huang 2009). Because the ensemble means of IMFs may not satisfy the condition of IMF, a post-processing EMD (Wu & Huang 2009) was applied on the decomposed components to guarantee that the final result satisfied the IMF criteria.

The light curve, x(t), was decomposed into 10 IMFs, denoted as c₁ to c₁₀, and a residual, r₁₀. The completeness of the IMFs is guaranteed, so that the light curve can be represented as

$$\begin{equation} x(t)=\sum _{j=1}^{10}c_j+r_{10}. \end{equation} \tag{ 1 }$$

The QPO signal lies within c₅, which can be examined by its variability timescale and its timing behavior such as its Lomb–Scargle periodogram and dynamic power spectrum. The orthogonality between neighboring components c_f and c_g can be expressed using the orthogonality index OI_fg (Huang et al. 1998):

$$\begin{equation} {\rm OI}_{fg}=\sum _{t}\frac{c_fc_g}{c_f^2+c_g^2}\protect. \protect \end{equation} \tag{ 2 }$$

The orthogonality index between c₄ and c₅ in this decomposition is OI₄₅ = 0.033, whereas the orthogonality index between c₅ and c₆ is OI₅₆ = 0.007. This means that the QPO signal is mostly concentrated within the IMF c₅, and the power barely leaks out to neighboring IMFs. Figure 2 shows the original light curve, EEMD high-pass-filtered light curve, QPO component, and long-term modulation. The EEMD high-pass-filtered light curve is a summation from c₁ to c₄, which represents the high-frequency noise in the data. The long-term modulation light curve, which represents the variability with timescale longer than the QPO, is a summation from c₆ to the residual that contains no extrema.

Zoom In Zoom Out Reset image size

We then applied the normalized Hilbert transform (Huang et al. 2009) to the QPO component to obtain the instantaneous frequency and amplitude. The instantaneous amplitude a_j(t) can be defined by the upper envelope of the absolute value of an IMF c_j(t). The Hilbert transform of a normalized IMF X_j(t) = c_j(t)/a_j(t) can be represented as

$$\begin{equation} Y(t)=\frac{1}{\pi }P\int _{-\infty }^{\infty }\frac{X(t^{\prime })}{t-t^{\prime }}dt^{\prime }\protect, \protect \end{equation} \tag{ 3 }$$

where P indicates the Cauchy principal value. Therefore, an analytical signal Z_j(t) and the angular phase function θ(t) can be defined as

$$\begin{equation} Z_j(t)=X_j(t)+iY_j(t)=e^{i\theta _j(t)}\protect. \protect \end{equation} \tag{ 4 }$$

Then, the instantaneous angular frequency ω_j(t) is defined as the time derivative of θ_j(t). As a result, the data can be expressed as

$$\begin{equation} x(t)=\sum _{j=1}^{10}a_j(t)\exp \left(i\int \omega _j(t)dt \right) +r_{10}\protect. \protect \end{equation} \tag{ 5 }$$

The variation of frequency and amplitude was plotted on a three-dimensional map, as shown in Figure 3, to demonstrate the drift of the QPO frequency and the variation of the QPO amplitude. We also produced the dynamic Lomb–Scargle periodogram for comparison. The window size of the dynamic power spectrum was set to 10 ks, and the moving step was 100 s. Instead of the original light curve, we applied the dynamic Lomb–Scargle algorithm on the detrended light curve, which was obtained by subtracting the long-term modulation light curve (Figure 2(d)) from the original light curve. Because the window size is approximately the same as the timescale of the flux-drop event ∼40–55 ks mentioned in Middleton et al. (2009), the power spectra are dominated by the long-term trend, and the QPO signals are strongly depressed when the window moves around here. The dynamic power spectrum is shown as the contour stacked on top of the Hilbert spectrum in Figure 3.

Zoom In Zoom Out Reset image size

From Figure 3, we found that the Hilbert spectrum agrees with the dynamic power spectrum. Both of them can describe the drift of the QPO frequency for timescales longer than ∼10 ks. The flux-drop event, together with a frequency increase in ∼40–55 ks, can be observed in this time–frequency map. Moreover, the Hilbert spectrum shows more details in both the time and frequency domains. The HHT provides not only the instantaneous frequency and amplitude but a well-defined phase function of QPO. We obtained a folded light curve by folding all the data according to the HHT phase. Figure 4 shows the light curves folded by the HHT phase, as well as ones folded by the cycle length and by a fixed period of 3733 s. The cycle length is the time difference between the neighboring fiducial points, which is defined by the local minima of the EEMD low-pass-filtered light curve obtained from a summation of QPO signal and the long-term modulation light curve. This is similar to the definition of local minima in Trowbridge et al. (2007), but we used the EEMD as a local filter instead of a Gaussian filter, and we used the local minima without fitting them with any functions. Because the QPO period is not stationary, the light curve folded with a fixed period of 3733 s would be expected to have distortion in its modulation shape. From Figure 4, the folded light curves produced by the HHT phase and the cycle length are similar, whereas the one folded by a fixed period is significantly different. We believe that the light curve folded by the HHT phase and that folded by the cycle length are better to describe the QPO profile than that folded by the fixed period.

Zoom In Zoom Out Reset image size

To study the evolution of the QPO period in detail, we applied the O − C analysis on the QPO cycle length. The variation of QPO period in RE J1034+396 is so dramatic that it can be detected in both the dynamic power spectrum and Hilbert spectrum. However, it is still worth tracing the evolution of the QPO period using O − C as long as the modulation is clearly seen without large observation gaps, which may cause the cycle count ambiguity. The fiducial point of each QPO cycle is defined as the one for the cycle length described in Section 3.1. This low-pass-filtered light curve provided a model-independent way to define the fiducial point, i.e., the minimum, in each cycle for the O − C analysis. We further applied a 10⁴ times Monte Carlo simulation to determine the error of the occurrence time of fiducial points. Then, we applied a linear ephemeris with the time of the first minimum as the phase zero epoch, and the computed period was 3733 s. Figure 5(a) shows the phase evolution of the QPO, and Figure 5(b) shows the corresponding light curve, EEMD low-pass-filtered light curve, and local minima. From the O − C results, we found that it is better to divide the evolution of the QPO into three epochs rather than two segments as proposed by Gierlínski et al. (2008). Epoch A is from the beginning of the observation to ∼42 ks, and the corresponding period obtained by fitting the phase evolution with a straight line is P_A = 4082.7 ± 24.6 s. Epoch B starts at ∼42 ks and ends at ∼55 ks, which corresponds to the flux-drop event. The period of epoch B obtained by O − C is P_B = 3094.4 ± 40.3 s, which is much shorter than that of epoch A. This decrease in period is also observed in the time–frequency map (see Figure 3). Because the flux-drop occurs with an increase in frequency, it strongly implies that this event is caused by the internal change in the QPO status rather than an external occultation as suggested by Middleton et al. (2009). Epoch C is t ≳ 55 ks, and the period is derived as P_C = 3784.1 ± 43.1 s, which is also significantly shorter than that of epoch A. The overall $\chi _{\nu }^2$ is 2.24, which implies a moderate fitting. Because the cycle length varies dramatically even between neighboring cycles, the large $\chi _{\nu }^2$ can be interpreted as resulting from the instability of the QPO period. We also fitted the segment 2 data with a straight line and found a period of P(t > 25 ks) = 3722.4 ± 27.1 s, which is consistent with the 3733 s determined by Gierlínski et al. (2008). However, the large $\chi _{\nu }^2$ value of 15.9 indicates an unacceptable fitting.

Zoom In Zoom Out Reset image size

The periodicities obtained by O − C can be further examined by the Lomb–Scargle periodogram. The power spectra of the three epochs are shown in Figures 5(c)–(e). The power spectrum of epoch A shows a strong peak at P_A = 4086.0 ± 10.6 s, whereas the one for epoch C also shows a strong peak at P_C = 3795.8 ± 16.2 s. The periods determined for epochs A and C are consistent with the result obtained by O − C. These periods are also consistent with the peaks in the power spectra of the entire data set (Figure 1). Furthermore, the periods of epochs A and C significantly deviate from each other at the ∼5σ and ∼10σ levels for the results obtained from the O − C analysis and Lomb–Scargle periodograms, respectively. The periodogram of epoch B exhibits a weak broad peak at P_B = 3282.2 ± 51.4 s. This period marginally agrees with that determined by O − C at only the 2σ level. This is acceptable because epoch B only consists of three cycles, and the QPO period changes between neighboring cycles.

Time lags either between different X-ray bands or between X-ray, UV, and optical wavelengths are intriguing issues of AGNs. Furthermore, time lags in the different timescales indicate different physical mechanisms. Middleton et al. (2011) found a soft lag in the folded QPO light curve. However, the QPO period was unstable, which implies that the phase lags are probably an artificial effect due to folding the non-stationary light curve with a fixed period, as we have shown in Section 3.1. We investigated this issue with the help of the HHT and O − C results.

We first reexamined the results obtained in Middleton et al. (2011). The light curve was divided into four energy bands: band 1 (0.2–0.3 keV), band 2 (0.3–0.5 keV), band 3 (0.5–1.0 keV), and band 4 (1.0–10.0 keV). We then folded the light curves of these four energy bands using a period of 3733 s and calculated the phase lags of bands 2–4 with respect to band 1 using the cross-correlation method. The folded light curves with the largest phase lags, i.e., bands 1 and 4, are shown in Figure 6(a), where the soft lag can be easily seen. Figure 6(b) shows the relationship between the phase lags and energies, which is similar to that obtained by Middleton et al. (2011). However, from both the Lomb–Scargle periodogram and time–frequency analysis in the previous sections, we noticed that a single period of 3733 s is insufficient to describe the QPO behavior. For example, we also folded the light curve using a 4083 s period, which is another significant peak in the Lomb–Scargle periodogram in the entire data set. The folded light curves of bands 1 and 4 are shown as Figure 6(c), and the corresponding phase lags are shown as Figure 6(d). No significant soft lag can be observed.

Zoom In Zoom Out Reset image size

Because the QPO period experienced three evolutionary epochs, we may obtain the phase lags between different energy bands of individual epochs. We first folded the light curves of four bands of epoch A using P_A = 4083 s, and then observed the phase lags between the folded light curves. No significant phase lags were detected between bands 1 and 2 or bands 1 and 3 (Figure 7(a)). The cross-correlation analysis shows that the folded light curve in band 4 slightly lags band 1, but the significance is only at ∼1σ level. We processed similar analysis for epoch C by folding the light curve with period of P_C = 3784 s. Figure 7(b) shows that the result is similar to that of epoch A, but band 4 appears to slightly lead band 1, although the significance is also at ∼1σ level. Finally, we folded the entire data set using the phase defined by the HHT to check the phase lags between different bands. However, no significant phase lags were detected (Figure 7(c)).

Zoom In Zoom Out Reset image size

We then examined the relationship between the flux and QPO period for all the data and three epochs, as defined in Section 3.2. To obtain an independent measurement of the period and to avoid the influence of intra-wave frequency modulation, i.e., the frequency variation within one cycle, we used the cycle length obtained in Section 3.2 as the period of each QPO cycle. On the other hand, we determined the flux of each QPO cycle by averaging the corresponding count rate of a long-term modulation light curve, which is shown in Figure 2(d), to decrease the influence of the noise and QPO modulation. Twenty QPO cycles were used in the correlation analysis. The relationship between the mean flux and QPO period is shown in Figure 8. The linear correlation coefficient of the entire data set is 0.67 with a null hypothesis probability of 3.4 × 10⁻⁴. This indicates a strong correlation between the mean flux and QPO period. For the individual epochs, epoch B contains only three cycles, so we ignore it and focus on the other two epochs. The linear correlation coefficient of epoch A is 0.11 with a null hypothesis probability of 0.68, which indicates no significant correlation. In contrast, the linear correlation coefficient of epoch C is 0.85 with a null hypothesis probability of 7.4 × 10⁻³, which is much stronger than that of epoch A.

Zoom In Zoom Out Reset image size

We further examined the slope between the QPO period and count rate to compare our results with that from Czerny et al. (2010). We made a logarithmic plot of the relationship between the mean count rate and QPO cycle length and then fitted the data points of the entire data set with a straight line. The slope is 0.30 ± 0.08, which is marginally consistent with the results of Czerny et al. (2010; 0.42 ± 0.05). We calculated the slopes for epochs A and C to be 0.09 ± 0.26 and 0.43 ± 0.12, respectively. The interpretation of the slopes and the dramatic change in the QPO period is discussed with the current QPO models in Section 5.

To determine the best spectral model, we divided the QPO phase into a high-intensity phase (0.2–0.55) and low-intensity phase (0.65–1.0) according to the folded light curve using the HHT phase. All the X-ray photons in the transition were discarded. Because the image of RE J1034+396 was highly concentrated and had a heavy pileup in the XMM-Newton PN and MOS chips (Middleton et al. 2009; Maitra & Miller 2010), excising the X-ray photons within the central region is needed. For the PN data of the entire data set, the 0.5–2 keV observed-to-model fraction of single events after excising the inner 32'' is 0.983 ± 0.017, which is acceptable at 1σ level (Maitra & Miller 2010). We also used the SAS task epatplot to examine this fraction of both the high- and low-intensity phases. The observed-to-model fraction of single events is 0.986 ± 0.014 after excising the inner 15'' region for the low-intensity state, whereas the fraction is 0.981 ± 0.021 after excising the inner 26'' region for the high-intensity state. To reduce the statistical error, we grouped the X-ray photons using grppha with a minimum of 100 photons in each energy bin. The XSPEC v12.8.0 of the HEASOFT package was used to perform the spectral fitting. The galactic hydrogen column density was fixed at 1.31 × 10²⁰ cm⁻²; however, an additional absorption to represent the absorption around the AGN was also set by adding a zphabs with a fixed redshift of 0.042. We then used the spectral model proposed by Middleton et al. (2009). A low-temperature Comptonized disk plus a hard power-law tail (hereafter known as compTT+pl) was used to fit both the high- and low-intensity phase data. Furthermore, we also added an additional blackbody component to the compTT+pl as the seed photons, hereafter known as bb+compTT+pl. Table 1 lists the best-fit parameters of these two models for both the high- and low-intensity phases.

Table 1. Best-fit Parameters for High- and Low-intensity Phases

Phase	High-intensity		Low-intensity
Model	compTT+pl	bb+compTT+pl	compTT+pl	bb+compTT+pl
kTBB (keV)	⋅⋅⋅	0.012 ± 0.007	⋅⋅⋅	0.013 ± 0.005
Γ	1.71 ± 1.1	2.60 ± 0.30	1.91 ± 0.52	2.22 ± 0.25
T0 (keV)	0.02 ± 0.08	=kTBB	0.02 ± 0.06	=kTBB
kTcomp (keV)	0.42 ± 0.59	0.16 ± 0.04	0.41 ± 0.28	0.24 ± 0.05
τ	7.1 ± 6.3	13.8 ± 3.3	7.05 ± 3.45	10.6 ± 1.6
$\chi _{\nu }^{2}$ (dof)	1.107(83)	1.124(82)	1.127(116)	0.982(115)

Notes. The kT_BB is the blackbody temperature of the blackbody model. For the Comptonized model, T₀ is the input seed photon temperature, kT_comp is the plasma temperature, and τ is the plasma optical depth. The power-law index of the power-law tail is Γ.

Download table as:  ASCIITypeset image

The fitted results reveal a blackbody emission of kT ∼ 10–15 eV, which probably represents the tail of the thermal emission from the gas in the inner accretion disk. Treating this blackbody emission as the seed photons for the bb+compTT+pl model, we obtained a good fit for both the high- and low-intensity phases with $\chi _{\nu }^2\approx 1.1$ and 1.0, respectively. The compTT+pl model, by comparison, cannot constrain a meaningful soft photon temperature in both the high- and low-intensity phases. Thus, we chose the bb+compTT+pl model to fit the phase-resolved spectra.

We further divided the X-ray photons into 16 phase bins according to the HHT phase. All the X-ray photons collected in each phase bin were then fitted by the best-fit model in Section 4.1. To increase the number of energy bins but keep sufficient statistics for each of the bins, especially the hard power-law tail with energy higher than 2 keV, we grouped the X-ray photons with a minimum of 25 photons in each energy bin. All the spectral parameters were first set to be free, and we found that they all have no significant variation during the QPO cycle except for the normalizations. Then, we froze the spectral parameters to the values obtained in Section 4.1, except for the normalizations, to estimate the flux variation of the spectral components during one QPO cycle. Because of the lack of sufficient hard X-ray photons with energies higher than 6 keV in individual phase bins, we calculated the flux between 0.3 and 6.0 keV. The fluxes of the blackbody component cannot be significantly detected because the normalization term of the blackbody component cannot be well constrained.

The left panel of Figure 9 shows the variation of those three spectral components. The power-law tail is more variable than the Comptonized component. We obtained the flux variation of these two components, shown in the right panel of Figure 9. We found that the fractional rms amplitude of the power-law component is 20.9%, which is approximately three times larger than that of the Comptonized component (6.6%). Furthermore, the phase lag between the Comptonized disk component and the power-law tail is −0.017 ± 0.030, where the error is estimated using a 10⁴ times Monte Carlo simulation, which indicates no significant phase lag. The lack of significant phase lag implies that if the Comptonized component is the reprocessed X-ray from the power-law component, the light-travel scale should be negligible.

Zoom In Zoom Out Reset image size

The most straightforward model to understand the QPO is that it is caused by the Keplerian motion of a temporary hot spot on the accretion disk (Pecháček et al. 2006; Bachetti et al. 2010). In this scenario, the QPO frequency is related to the radius of Keplerian orbit. From the HHT and O − C analysis, we found that the QPO period changes dramatically even between neighboring cycles. This indicates that the hot spot is not strictly anchored on the accretion disk. Instead, the hot spot is wobbling around a characteristic radius. Furthermore, the evolution of the QPO period can be described by three epochs in this observation. This evolution can be interpreted as the change of characteristic radius caused by the instability of the accretion disk.

From the analysis of the relationship between modulation period and the flux, we observed a positive correlation between the QPO period and the flux. This relationship, which basically agrees with that observed in Czerny et al. (2010), does not agree with the prediction of this model. However, some indications of this model are still worth addressing. Czerny et al. (2010) compared the fractional rms amplitude of RE J1034+396 and that of the light curve calculated by a spotted accretion disk (Pecháček et al. 2006). They concluded that the inclination angle should be 2°–3° if the hot spot is close to the central black hole. However, not only the flare but also the disk, corona, and probably other materials near the black hole contribute to the X-ray emissions. For example, we obtained three components in the spectral fitting. Those X-ray sources would contribute non-modulated or lower-amplitude X-ray emissions such that the fractional amplitude and the inclination angle would be underestimated.

Besides the fractional amplitude, another possible indicator that can be used to estimate the inclination angle is the modulation shape. From the simulated light curve in Czerny et al. (2010) and Pecháček et al. (2006), the QPO modulation shape became highly non-sinusoidal for higher inclination angles. A non-sinusoidal light curve can be expressed as harmonics in the Fourier analysis; however, we could not observe any significant harmonics in the power spectrum. There are two possible explanations for this. The first is that the inclination angle of RE J1034+396 is low, and the QPO profile is close to sinusoidal such that the harmonics are not significantly detected. Another possibility is that the QPO period varies so dramatically that the harmonics are more strongly depressed than the fundamental. Fortunately, the HHT provides an indicator to probe the non-sinusoidal level of the data. Unlike harmonics in a Fourier analysis, the non-sinusoidal light curve would produce intra-wave frequency modulations, i.e., the frequency variation within one cycle (Huang et al. 1998) in the instantaneous frequency. The frequency derived from the cycle length in the O − C analysis can be treated as the inter-wave frequency modulation, i.e., the frequency variation between cycles. Thus, the intra-wave frequency modulation can be estimated by taking the difference between the instantaneous frequency and inter-wave frequency modulation. We then created simulated light curves according to Equations (4) and (10) in Pecháček et al. (2006) for different black hole masses and inclination angles. The orbital radius of the hot spot is calculated using Keplers law for orbital periods of 3784 and 4083 s. The amplitude of the intra-wave frequency modulation of a simulated light curve can also be calculated using the Hilbert transform on the simulated light curve. For a mass range (1–4) × 10⁶ M_☉, we can estimate the inclination angle by matching the intra-wave frequency modulation amplitude of the simulated light curves with the observed one. Figure 10 shows the relationship between the estimated inclination angle and black hole mass. We found that the inclination angle was between ∼27° and ∼57°. This estimation of inclination angle is probably an overestimate because the noise in the light curve would also contribute to the intra-wave modulation.

Zoom In Zoom Out Reset image size

Relativistic diskoseismology, i.e., the oscillations that are trapped in the inner accretion disk, may also be responsible for the presence of QPOs (Perez et al. 1997; Silbergleit et al. 2001; Wagoner et al. 2001; Ortega-Rodríguez et al. 2002). Three oscillation modes, i.e., g, c, and p modes, are investigated. The fundamental frequency of these three modes mainly depends on the mass of the black hole, i.e., f∝1/M_BH, where f is the QPO frequency and M_BH is the mass of the black hole. From Wagoner et al. (2001), the g mode is dominated by the gravitational–centrifugal force, in which the oscillation frequency anti-correlates with the flux of the accretion disk. The c mode is caused by the non-radial corrugation on the inner accretion disk, while the oscillation frequency is highly dependent on the spin of the black hole. The p mode is dominated by the restoring force due to the pressure gradient, in which the oscillation frequency shows a positive correlation with the flux of the disk. Only the p mode shows a positive correlation between the flux and QPO frequency, but it is the opposite of our results. Thus, the QPO in RE J1034+396 is not likely caused by the p-mode oscillation. For a black hole with mass of 10 M_☉, the g-mode oscillation is between ∼70 and ∼110 Hz for different angular momenta and luminosities (see Figure 1 of Wagoner et al. 2001). If we enlarge the black hole mass to that in RE J1034+396, the g-mode oscillation period is between ∼900 and ∼1400 s for a 1 × 10⁶ M_☉ black hole, or between ∼3600 and ∼5700 s for a 4 × 10⁶ M_☉ black hole. Thus, the preferred black hole mass is close to the upper limit of the estimation made by Bian & Huang (2010) if we assume that the QPO is caused by the g-mode oscillation. However, it is difficult for the slope between the g-mode oscillation period and luminosity to exceed ∼0.1 unless the angular momentum of the black hole is extremely high (see Equation (3) in Wagoner et al. 2001). Thus, the g-mode oscillation cannot provide the observed slope (∼0.3), although this value is less reliable relative to the correlation coefficient because the observed X-ray flux may also contain the emissions from other components. On the other hand, the c-mode oscillation can produce a larger slope between the flux and QPO frequency in the range of our observed value, but the angular momentum of the black hole should be close to zero (see Equation (6) in Wagoner et al. 2001). However, the fundamental c-mode frequency of such a slow-rotating black hole is very low. For example, a 10 M_☉ black hole with angular momentum of a = 0.05 has a fundamental c-mode oscillation frequency of f ≲ 1 Hz. This oscillation corresponds to a period longer than 10⁵ s for a 10⁶ M_☉ black hole. Thus, if the c-mode oscillation is responsible for the observed QPO, it should be high-order harmonics.

We also studied the relationship between the QPO frequency and flux in detail by dividing the evolution of the QPO into three epochs. A change in correlation between epochs was observed: a significant anti-correlation between the frequency and flux was observed in epoch C, but no significant correlation was found in epoch A. This is possibly caused by a change in oscillation mode. A correlation in epoch A was not significantly detected; however, we found that the first three cycles appear to behave similar to the flux-drop event in epoch B. If we ignore these three cycles, the remaining cycles in epoch A show a strong anti-correlation between the cycle length and flux with a linear correlation coefficient of −0.61. Thus, the QPO in epoch A is probably caused by the p-mode oscillation. After epoch B, i.e., the flux-drop event, the g-mode oscillation dominated the observed QPO phenomenon such that we observed a significant direct correlation between the cycle length and flux.

All the discussions above are based on the assumption that the luminosity is less than the Eddington luminosity. This is a necessary condition for the thin-disk assumption. However, the luminosity of RE J1034+396 is probably super-Eddington (Gierlínski et al. 2008; Middleton et al. 2009). Thus, we could not exclude all the possible oscillation modes in a super-Eddington accretion disk.

Das & Czerny (2011) proposed a novel model of QPO to investigate the observed slope between the QPO period and X-ray flux in RE J1034+396. In their model, QPO is caused by the oscillations of the shock formed inside the hot accretion flow. The variation in the QPO period is interpreted as the drifting of the shock location. This model successfully interpreted the correlation between the QPO period and X-ray flux in Czerny et al. (2010) with slope of 0.92 ± 0.03. However, we have demonstrated that the slope changes as QPO evolves in this observation. In fact, even the relationship likely changes. For example, the QPO period and flux exhibit no significant correlation in epoch A, or even anti-correlation when we omit the first three cycles that show another possible flux-drop event. Thus, the entire observation cannot be interpreted by this single scenario.