Phase-resolved Analyses of Millihertz Quasi-periodic Oscillations in 4U 1636-53 using the Hilbert–Huang Transform

Various kinds of quasi-periodic oscillations (QPOs) that exhibit boxer peaks in the power spectra have been observed in many low-mass X-ray binaries (LMXBs). These QPOs are believed to reflect phenomena associated with accretion disks (see van der Klis 2006 for an extensive review). However, a new kind of QPO, referred to as a millihertz (mHz) QPO, was discovered in the accreting neutron-star binaries 4U 1636-53 and 4U 1608-52 and has occasionally been observed in Aql X-1 (Revnivtsev et al. 2001). Revnivtsev et al. (2001) suggested that these mHz QPOs are caused by a special mode of nuclear burning on the neutron-star surface because the flux in 4U 1608-52 was in a narrow range during mHz QPOs, which is approximately consistent with the flux in which type-I X-ray bursts cease to exist. However, they did not exclude the possibility that the mHz QPOs were caused by instability of the accretion disk. Revnivtsev et al. (2001) also discovered that all mHz QPOs occur before a type-I burst and disappear afterward. Following this discovery, mHz QPOs were also detected in numerous LMXBs, such as 4U 1323-62 (Strohmayer & Altamirano 2012), IGR J17480-2446 in Terzan 5 (Linares et al. 2012), and "clocked" burster GS 1826-238 (Strohmayer et al. 2018). Yu & van der Klis (2002) suggested that the kilohertz (kHz) QPO frequency in 4U 1608-52 is anticorrelated with the count rate during a mHz QPO cycle. Because the kHz QPO is related to the Keplerian orbital frequency at the inner edge of the disk, a positive correlation between mHz QPO and kHz QPO should have been detected if the mHz QPO was produced from the modulation of the disk. This finding provides evidence that mHz QPOs comprise a phenomenon that occurs on the neutron-star surface instead of in the accretion disk. Altamirano et al. (2008) found that the QPO frequency decreases with time during a mHz oscillation until it disappears as a type-I X-ray burst occurs in 4U 1636-53.

Previous studies showed that mHz QPOs are associated with nuclear burning on a neutron-star surface (Revnivtsev et al. 2001; Altamirano et al. 2008). The fuel for nuclear reaction is obtained from accretion. Therefore, the accretion rate is an important factor for nuclear burning on the neutron-star surface. Stable burning indicates that the nuclear burning rate (ε_nuc) and cooling rate (ε_cool) are balanced, and both depend on temperature and the depth of the fuel layer (Fujimoto et al. 1981; Keek & Heger 2015). Fuel piled up on the layer is compressed by a strong gravitational force, which increases the temperature and density in the fuel layer. When the ignition condition (ε_nuc > ε_cool) is reached, a thermonuclear runaway would occur. Thermonuclear runaway causes the luminosity to increase in a short time. This phenomenon is called a type-I X-ray burst. Most sources of type-I X-ray bursts have mixed fuel (hydrogen and helium) on the neutron-star surface, and helium burning is ignited by hydrogen burning. For a high accretion rate ($\gt 11 \% {\dot{M}}_{\mathrm{Edd}}$), the accretion provides sufficient fuel to the burning layer. Hydrogen burning dominates in the fuel layer, and helium accumulates at the bottom of the layer. As helium burning is ignited, the 3α reaction releases more energy and creates more carbon, which enhances the hot CNO cycle. Hydrogen burning also promotes the 3α reaction. This burning layer increases ε_nuc further, and then runaway thermonuclear reactions occur while reaching the ignition condition (Galloway & Keek 2017).

Heger et al. (2007) proposed that the mHz QPO occurs in a special condition between stable and unstable burning, which is referred to as marginal burning. The nuclear reaction rate increases as the temperature in the fuel layer increases, which causes the fuel in the layer to be consumed faster than it is supplied by accretion. Thus, the thickness of the burning layer, y, decreases. The cooling rate is proportional to F/y, where F is the outward flux. The outward flux F ≈ acT⁴/3κy, where a is the radiation constant, c is the speed of light, and κ is the opacity (Bildsten 1998). Therefore, as the thickness of the burning layer y decreases, the cooling rate increases to reduce the temperature and the nuclear reaction rate. As the rate at which fuel is consumed by the nuclear reaction is less than the rate at which it is supplied by accretion, the thickness of the burning layer y increases, which leads to a lower cooling rate. Therefore, the temperature and nuclear reaction rate increase again. Heger et al. (2007) estimated that it takes about $P\cong \sqrt{{t}_{\mathrm{therm}}{t}_{\mathrm{acc}}}\sim 120\,{\rm{s}}$ to complete this QPO cycle. According to this model, marginally stable nuclear burning occurs at an accretion rate close to the Eddington rate (Bildsten 1998; Heger et al. 2007). However, the results of observations showed that the accretion rates are much lower than the Eddington rate in 4U 1608-52 and 4U 1636-53 (Revnivtsev et al. 2001). Keek et al. (2009) attempted to find a possible solution in mixing processes. Furthermore, the width of the interval in accretion rate between stable and unstable nuclear burning depends strongly on the composition of the layer and reaction rates (Keek et al. 2014).

Lyu et al. (2014, 2015) analyzed the mHz QPOs in 4U 1636-53 using XMM-Newton observations. They also observed the drift in the QPO frequency, which is positively correlated with the temperature of the blackbody component from the neutron-star surface. Such drifts could be caused by cooling of the burning layers. By analyzing 39 type-I X-ray bursts and the mHz QPO in 4U 1636-53, Lyu et al. (2016) determined that the mHz QPO can only be detected when the convexity of the corresponding type-I X-ray burst is positive. According to the model proposed by Maurer & Watts (2008) and Cooper & Narayan (2007), the convexity links to the ignition site of the bursts on the neutron-star surface. Lyu et al. (2016) suggested that the marginally stable nuclear burning process may occur at the neutron-star equator. Through phase-resolved spectroscopy of mHz oscillations of 4U 1636-53, Stiele et al. (2016) concluded that the mHz QPOs are caused by a variable surface area of nuclear burning on neutron-star surface with constant temperature. Nevertheless, Strohmayer et al. (2018) detected the oscillation of the blackbody temperature component in an effectively constant emitting area in the 8 mHz QPO of GS 1826-238.

Phase-resolved analysis may help us to further understand the nature of mHz QPOs. However, the definition of the phase with traditional analysis methods (e.g., Fourier transform) for unstable oscillation phenomena, such as QPOs, is difficult. In this research, we attempt to use Hilbert–Huang transform (HHT) to analyze mHz QPOs. The HHT is a powerful tool for analyzing phenomena with nonstationary periodicity and has been successfully applied in astronomical research, such as for the QPO in the active galactic nucleus RE J1034 + 396 (Hu et al. 2014) and the 4 Hz QPO in the black hole X-ray binary XTE J1550-564 (Su et al. 2015). The HHT enables us to not only trace the variation in frequency in the QPO but also to process phase-resolved analyses even though the periodicity is unstable. In this paper, we present our analysis of the evolution of the frequency of mHz QPOs in 4U 1636-53 and the QPO phase-resolved spectral variations. In Section 2, we briefly introduce the XMM-Newton observations and data reduction process. In Section 3, we describe how HHT analysis is applied in mHz QPOs and demonstrate the phase-resolved spectral analysis that reveals the variations in spectral parameter for the complete cycle. From the correlation between spectral parameters and fluxes, we find that the oscillation is mainly attributed to variations in the area emitting blackbody radiation in three out of four observations with mHz QPO detections, whereas the other one shows a concurrent variation of temperature and flux with a constant emitting area. Finally, we discuss the HHT-based timing property and possible implications for our spectral analysis in Section 4.

The data being analyzed with the HHT in this work were collected by XMM-Newton between 2000 and 2015, a total of 12 observations as listed in Table 1. The HHT favors the data detected by high-altitude X-ray telescopes, such as XMM-Newton because the gaps in observation from low-altitude satellites would induce too much aliasing. For 4U 1636-53, XMM-Newton can provide ∼30,000 s continuous exposures, which are sufficient for the HHT analysis of ∼8 mHz QPOs. The data analyzed in this work were collected by the European Photon Imaging Camera (EPIC-PN) of XMM-Newton, and the data reduction was performed by the Science Analysis System (SAS) version 16.0.0. The current calibration file (CCF) of XMM-Newton in this work is updated to XMM-CCF-REL-371. Similar to the selection criteria proposed by Lyu et al. (2015), we selected a region 10 columns wide that is centered at the position of the source and only single and double events (PATTERN ≤ 4). The xmmselect, the graphical interface to the SAS extractor, produced event files comprising photon energies of 0.8 to 11 keV with a time resolution of 0.03 ms in the timing mode.

3.1. Hilbert–Huang Transform Analysis

All photon arrival times were first corrected to the barycenter of the solar system by the SAS task barycen. In order to optimize for the HHT analysis of the ∼8 mHz QPO, we binned the events collected in all 12 XMM-Newton observations into light curves with a resolution of 1.3333 s. The selection of this time resolution will be explained later in this section. The dynamic power spectral analysis was adopted for these 12 light curves to determine the time intervals with significant ∼8 mHz QPO detections. These power spectra were obtained by using the Lomb–Scargle power periodogram (Scargle 1982) with a window size of 1130 s, ∼10 cycles of QPO variations, and a moving step size of 100 s. The ∼8 mHz QPO can be clearly seen only in 4 of 12 light curves, labeled as Obs1 to Obs4 (see Table 1), whose dynamic power spectra are shown in Figure 1. Except for Obs3, the frequency was decreased before a type-I X-ray burst according to the dynamic power spectrum, consistent with the previous study (Lyu et al. 2015).

Zoom In Zoom Out Reset image size

Although the variation in the frequency and power can be clearly seen, the resolution of the dynamic power spectrum is limited by the window size, that is, δfδT ≈ 1, where δf and δT are the resolutions in frequency and time (i.e., window size), respectively. This limitation prevents us from further analysis of this variation. To investigate the ∼8 mHz QPO in more detail, the HHT was applied in the rest of this work. However, a continuous light curve is preferred for the HHT analysis but there are gaps in Obs2 and Obs3. Therefore, we only analyzed the parts after the gaps because they are closer to the corresponding type-I bursts.

We performed HHT analysis for these four selected light curves. The HHT is a method for analyzing nonlinear and nonstationary signals. The method consists of two major steps (Huang & Wu 2008): (1) using empirical mode decomposition (EMD) to adaptively decompose a time series into intrinsic mode functions (IMFs). (2) Obtaining the instantaneous frequency and amplitude of the IMFs through the Hilbert transform. EMD is a method for decomposing any time series into several IMFs. This process was developed by Huang et al. (1998) and then further improved by Huang et al. (2009). IMFs are time series that are suitable for the Hilbert transform to resolve the instantaneous frequency and amplitude of nonstationary periodic oscillations, such as QPOs.

In this work, we applied the fast-complementary ensemble empirical mode decomposition (CEEMD) method with post-processing (Huang et al. 2009; Yeh et al. 2010; Wang et al. 2014), a modified version of EMD. EMD is a sifting process for separating oscillation modes from original data by subtracting the local means in the data (Huang et al. 1998). These decomposed components are IMFs that satisfy the following conditions: (1) the number of extrema and the number of zero crossings must either be equal or differ by at most one, and (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero (Huang et al. 1998). However, IMFs extracted by EMD may suffer from the mode mixing problem, in which a modulation with the same timescale is distributed across different IMFs (Yeh et al. 2010). CEEMD can overcome this mode mixing problem by (1) adding white noise to original data, (2) decomposing the noisy data into IMFs by EMD, (3) repeating steps (1) and (2) several times with different white noise each time, (4) averaging these IMFs. Because CEEMD involves numerous summations of IMFs, it indicates that the CEEMD components of data may not be IMFs. In order to guarantee that the final result satisfied the IMF criteria, the post-processing EMD (Huang et al. 2009) was applied to the decomposed components.

The light curve, x(t), may be expressed as the sum of the IMFs, c_j(t), and the residual r_n(t), $x(t)={\sum }_{j=1}^{n}{c}_{j}(t)+{r}_{n}(t)$, where n is the number of IMFs. Because the variation in the timescale of the first IMF component is approximately three data points, i.e., 3 × 1.3333 s = 4 s, and the variation in the timescale of the second component is twice that of the first (Huang et al. 1998) and so on, the timescale in which we are interested (∼125 s) is expected to be concentrated at the sixth IMF component (4 × 2^(6 − 1) = 128 s). After decomposing the light curve, we determined that the ∼8 mHz oscillation can be observed in the sixth IMF, that is, c₆(t) as expected (see Figure 2), by confirming the orthogonality of the IMF components (Huang et al. 1998; Hu et al. 2014). We then applied the normalized Hilbert transform (Huang et al. 2009) to obtain the instantaneous phase, frequency, and amplitude of IMF c₆(t). The instantaneous amplitude a₆(t) is defined as the cubic Hermite spline envelope of the local maxima of the absolute values of the IMF c₆(t) (Huang et al. 2009). The Hilbert transformation of the normalized IMF X₆(t) = c₆(t)/a₆(t) can be represented as

$$\begin{eqnarray}&&{Y}_{6}(t)=\displaystyle \frac{1}{\pi }P{\int }_{-\infty }^{\infty }\displaystyle \frac{{X}_{6}(t^{\prime} )}{t-t^{\prime} }{dt}^{\prime} ,\end{eqnarray} \tag{ 1 }$$

where P is the Cauchy principal value. Thus, we can define an analytical signal Z₆(t) and the instantaneous phase function θ₆(t) as

$$\begin{eqnarray}&&{Z}_{6}(t)={X}_{6}(t)+{{iY}}_{6}(t)={e}^{i{\theta }_{6}(t)}.\end{eqnarray} \tag{ 2 }$$

Therefore, the instantaneous frequency ν₆(t) can be defined as

$$\begin{eqnarray}&&{\nu }_{6}(t)=\displaystyle \frac{1}{2\pi }\displaystyle \frac{d{\theta }_{6}(t)}{{dt}}.\end{eqnarray} \tag{ 3 }$$

Figure 3 shows Hilbert spectra of the ∼8 mHz oscillation with more detailed variations in both frequency and amplitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To significantly improve detection, the confidence limits can be determined by repeating the CEEMD 10³ times with different sets of white noise to generate a thousand versions of IMF c₆(t), and then the means and standard deviations of the instantaneous amplitude and frequency from these IMF sets can be calculated. Finally, a significant signal is defined as one with an amplitude above the average of the 3σ lower limit (∼6 counts s⁻¹) (see Su et al. 2015 for details); the corresponding data will be selected for further analysis.

3.2. Oscillation Profile and Phase-resolved Spectra

After the instantaneous QPO phases $\phi (t)=\mathrm{frac}[{\theta }_{6}(t)]$ are evaluated by the HHT (Equation (2)), the oscillation profile can be constructed by folding them even though the oscillation period is unstable. The oscillation profiles were constructed by binning the phases into 20 bins per cycle, as shown in Figure 4. The nonsinusoidal nature can be clearly seen in the oscillation profile, which may be approximately described as a fast rise of ∼0.3 cycle and an exponential decay of ∼0.5 cycle.

Zoom In Zoom Out Reset image size

With the well-defined phase (Equation (2)), phase-resolved spectral analysis can be performed. We employed evselect, a standard SAS tool, to extract energy spectra. The spectra were extracted according to the selection criteria that PATTERN ≤ 4 and FLAG = 0. Because the source is very bright, the whole CCD was contaminated by source photons (Lyu et al. 2014, 2015). The background spectra were evaluated by the black hole candidate GX339-4 in quiescent state, which is in a region of sky with a similar column density along the line of sight (Hiemstra et al. 2011; Sanna et al. 2013).

To obtain the phase-resolved spectra, we divided an oscillation cycle into 20 bins according to the phase defined by the HHT, as we did for the oscillation profile, and extracted their individual energy spectra. Using XSPEC v12.10.1, these spectra were fitted with the model (Lyu et al. 2015; Stiele et al. 2016):

$$\begin{eqnarray*}&&{\rm{PHABS}}\times ({\rm{BBODYRAD}}+{\rm{DISKBB}}+{\rm{NTHCOMP}}),\end{eqnarray*}$$

where PHABS is photoelectric absorption, BBODYRAD is blackbody emission from the neutron-star surface, DISKBB (Mitsuda et al. 1984; Makishima et al. 1986) is multicolor disk blackbody emission to describe the thermal emissions from an accretion disk at various temperatures, and NTHCOMP (Zdziarski et al. 1996; Życki et al. 1999) is the inverse Compton scattering process in the corona, with the seed thermal photons coming from the accretion disk (Sanna et al. 2013; Lyu et al. 2014). The interstellar absorption was fixed at 0.36 × 10²² cm⁻² (Sanna et al. 2013), selecting the solar abundance from Wilms et al. (2000) and the cross section from Verner et al. (1996). Because the mHz QPO is considered as metastable nuclear burning on the neutron-star surface, only the parameters in BBODYRAD were allowed to vary for different phase bins, with the parameters of other components being tied to the optimal values evaluated from the detected spectra of the whole cycle. All 20 phase bin spectra of the four observations are well fitted by this model, and therefore the variations in spectral parameter on the neutron-star surface can be resolved for the complete cycle. Figure 5 shows the typical spectra from the first phase bin of four observations.

Zoom In Zoom Out Reset image size

Figure 4 illustrates modulations of the two untied parameters in BBODYRAD, the temperature T_bb and the apparent area ${R}_{\infty }^{2}$, as well as the oscillation profile. The apparent area $({R}_{\infty }^{2})$ with the mHz QPO phase was calculated using the relation (Sztajno et al. 1985; Stiele et al. 2016)

$$\begin{eqnarray*}&&{R}_{\mathrm{BB}}^{2}={R}_{\infty }^{2}\times {f}_{\mathrm{col}}^{4}\times \left(1-2\displaystyle \frac{M}{{R}_{\mathrm{NS}}}\right),\end{eqnarray*}$$

where G = c = 1. The source surface (${R}_{\mathrm{BB}}^{2}$) was derived from BBODYRAD components with a distance of 6 kpc (Galloway et al. 2006). The factor f_col was set to 1.6 with a helium-enriched environment (Lyu et al. 2015). According to the bolometric flux oscillations that occur during the rise of the X-ray burst from Rossi X-ray Timing Explorer (RXTE) observations, M/R_NS is 0.126 (Nath et al. 2002). Tables 2 and 3 list the untied spectral parameters from BBODYRAD and tied parameters from DISKBB and NTHCOMP, respectively.

Table 2.  All BBODYRAD Parameters in This Work

Note	Phase	Temperature (keV)	Area Normalized	χ2/dof	Note	Phase	Temperature (keV)	Area Normalized	χ2/dof
Obs1	1	${0.5477}_{-0.0080}^{+0.0080}$	${789}_{-34}^{+35}$	$138.46/155$	Obs2	1	${0.5952}_{-0.0053}^{+0.0038}$	${665}_{-22}^{+19}$	$173.10/169$
Obs1	2	${0.5507}_{-0.0083}^{+0.0083}$	${757}_{-34}^{+34}$	$155.41/159$	Obs2	2	${0.5971}_{-0.0052}^{+0.0038}$	${639}_{-21}^{+19}$	$165.19/170$
Obs1	3	${0.5472}_{-0.0084}^{+0.0085}$	${764}_{-34}^{+35}$	$128.89/156$	Obs2	3	${0.5913}_{-0.0055}^{+0.0039}$	${644}_{-22}^{+19}$	$185.82/171$
Obs1	4	${0.5430}_{-0.0090}^{+0.0090}$	${767}_{-18}^{+36}$	$155.12/156$	Obs2	4	${0.5895}_{-0.0056}^{+0.0040}$	${647}_{-22}^{+19}$	$154.66/171$
Obs1	5	${0.5412}_{-0.0091}^{+0.0091}$	${756}_{-29}^{+36}$	$172.28/157$	Obs2	5	${0.5909}_{-0.0057}^{+0.0040}$	${625}_{-22}^{+19}$	$210.04/170$
Obs1	6	${0.5390}_{-0.0092}^{+0.0092}$	${761}_{-35}^{+36}$	$119.45/159$	Obs2	6	${0.5948}_{-0.0052}^{+0.0040}$	${606}_{-25}^{+19}$	$199.24/170$
Obs1	7	${0.5377}_{-0.0095}^{+0.0094}$	${762}_{-36}^{+36}$	$149.83/159$	Obs2	7	${0.5904}_{-0.0057}^{+0.0040}$	${618}_{-26}^{+19}$	$160.11/171$
Obs1	8	${0.5387}_{-0.0093}^{+0.0092}$	${748}_{-35}^{+35}$	$148.07/157$	Obs2	8	${0.5917}_{-0.0058}^{+0.0040}$	${604}_{-26}^{+19}$	$121.97/171$
Obs1	9	${0.5381}_{-0.0095}^{+0.0095}$	${754}_{-36}^{+36}$	$134.25/159$	Obs2	9	${0.5960}_{-0.0056}^{+0.0040}$	${594}_{-21}^{+19}$	$166.33/171$
Obs1	10	${0.5468}_{-0.0091}^{+0.0091}$	${709}_{-33}^{+34}$	$157.38/158$	Obs2	10	${0.5932}_{-0.0058}^{+0.0040}$	${599}_{-19}^{+19}$	$166.44/177$
Obs1	11	${0.5458}_{-0.0089}^{+0.0089}$	${721}_{-34}^{+34}$	$178.75/157$	Obs2	11	${0.5893}_{-0.0059}^{+0.0041}$	${606}_{-21}^{+19}$	$157.21/172$
Obs1	12	${0.5418}_{-0.0093}^{+0.0093}$	${737}_{-18}^{+35}$	$143.77/158$	Obs2	12	${0.5896}_{-0.0058}^{+0.0040}$	${606}_{-26}^{+19}$	$185.17/172$
Obs1	13	${0.5323}_{-0.0099}^{+0.0099}$	${774}_{-37}^{+38}$	$141.56/157$	Obs2	13	${0.5897}_{-0.0058}^{+0.0040}$	${612}_{-21}^{+19}$	$154.24/172$
Obs1	14	${0.5312}_{-0.0101}^{+0.0100}$	${793}_{-37}^{+39}$	$170.59/157$	Obs2	14	${0.5903}_{-0.0059}^{+0.0041}$	${607}_{-21}^{+19}$	$161.58/171$
Obs1	15	${0.5363}_{-0.0093}^{+0.0093}$	${769}_{-36}^{+36}$	$138.23/156$	Obs2	15	${0.5965}_{-0.0051}^{+0.0039}$	${592}_{-26}^{+19}$	$178.58/172$
Obs1	16	${0.5420}_{-0.0091}^{+0.0090}$	${741}_{-32}^{+35}$	$135.26/158$	Obs2	16	${0.5967}_{-0.0056}^{+0.0040}$	${596}_{-26}^{+19}$	$181.53/170$
Obs1	17	${0.5538}_{-0.0082}^{+0.0082}$	${731}_{-33}^{+30}$	$154.51/157$	Obs2	17	${0.5937}_{-0.0050}^{+0.0039}$	${633}_{-27}^{+19}$	$196.67/170$
Obs1	18	${0.5516}_{-0.0081}^{+0.0082}$	${765}_{-33}^{+34}$	$171.72/156$	Obs2	18	${0.5914}_{-0.0053}^{+0.0041}$	${662}_{-22}^{+19}$	$184.78/171$
Obs1	19	${0.5565}_{-0.0081}^{+0.0080}$	${748}_{-33}^{+30}$	$149.87/157$	Obs2	19	${0.5948}_{-0.0048}^{+0.0037}$	${658}_{-25}^{+19}$	$192.30/171$
Obs1	20	${0.5591}_{-0.0075}^{+0.0076}$	${747}_{-32}^{+32}$	$160.29/156$	Obs2	20	${0.5975}_{-0.0050}^{+0.0036}$	${661}_{-22}^{+19}$	$202.71/171$
Obs3	1	${0.5398}_{-0.0047}^{+0.0048}$	${714}_{-30}^{+30}$	$180.53/171$	Obs4	1	${0.5637}_{-0.0056}^{+0.0058}$	${820}_{-20}^{+23}$	$154.50/169$
Obs3	2	${0.5388}_{-0.0048}^{+0.0049}$	${697}_{-31}^{+30}$	$189.44/170$	Obs4	2	${0.5591}_{-0.0059}^{+0.0061}$	${818}_{-23}^{+23}$	$190.83/169$
Obs3	3	${0.5270}_{-0.0051}^{+0.0052}$	${737}_{-31}^{+31}$	$224.64/170$	Obs4	3	${0.5549}_{-0.0062}^{+0.0063}$	${825}_{-24}^{+24}$	$177.72/169$
Obs3	4	${0.5240}_{-0.0053}^{+0.0055}$	${737}_{-31}^{+31}$	$188.45/171$	Obs4	4	${0.5597}_{-0.0062}^{+0.0063}$	${776}_{-23}^{+23}$	$184.26/170$
Obs3	5	${0.5292}_{-0.0052}^{+0.0054}$	${701}_{-31}^{+30}$	$155.28/171$	Obs4	5	${0.5584}_{-0.0064}^{+0.0066}$	${767}_{-23}^{+23}$	$156.55/169$
Obs3	6	${0.5289}_{-0.0055}^{+0.0056}$	${676}_{-31}^{+30}$	$180.79/170$	Obs4	6	${0.5530}_{-0.0033}^{+0.0067}$	${796}_{-23}^{+24}$	$193.24/170$
Obs3	7	${0.5291}_{-0.0055}^{+0.0057}$	${668}_{-31}^{+30}$	$150.00/170$	Obs4	7	${0.5549}_{-0.0064}^{+0.0066}$	${777}_{-23}^{+23}$	$201.52/170$
Obs3	8	${0.5293}_{-0.0055}^{+0.0057}$	${659}_{-31}^{+30}$	$199.10/170$	Obs4	8	${0.5572}_{-0.0032}^{+0.0066}$	${767}_{-23}^{+23}$	$123.72/170$
Obs3	9	${0.5288}_{-0.0054}^{+0.0056}$	${677}_{-31}^{+30}$	$189.07/171$	Obs4	9	${0.5615}_{-0.0064}^{+0.0066}$	${734}_{-22}^{+23}$	$161.16/170$
Obs3	10	${0.5317}_{-0.0055}^{+0.0057}$	${652}_{-30}^{+30}$	$146.95/171$	Obs4	10	${0.5559}_{-0.0064}^{+0.0067}$	${763}_{-23}^{+23}$	$153.28/168$
Obs3	11	${0.5287}_{-0.0058}^{+0.0059}$	${661}_{-31}^{+30}$	$180.95/171$	Obs4	11	${0.5595}_{-0.0064}^{+0.0066}$	${740}_{-23}^{+23}$	$170.36/170$
Obs3	12	${0.5271}_{-0.0055}^{+0.0056}$	${672}_{-31}^{+30}$	$176.76/170$	Obs4	12	${0.5549}_{-0.0067}^{+0.0069}$	${754}_{-23}^{+24}$	$207.85/168$
Obs3	13	${0.5287}_{-0.0056}^{+0.0058}$	${660}_{-31}^{+30}$	$189.05/171$	Obs4	13	${0.5605}_{-0.0064}^{+0.0065}$	${737}_{-23}^{+23}$	$165.73/170$
Obs3	14	${0.5268}_{-0.0056}^{+0.0058}$	${683}_{-31}^{+31}$	$208.90/172$	Obs4	14	${0.5535}_{-0.0066}^{+0.0068}$	${767}_{-23}^{+24}$	$174.87/172$
Obs3	15	${0.5368}_{-0.0055}^{+0.0056}$	${635}_{-31}^{+30}$	$179.65/171$	Obs4	15	${0.5587}_{-0.0064}^{+0.0066}$	${755}_{-23}^{+23}$	$168.19/168$
Obs3	16	${0.5322}_{-0.0054}^{+0.0055}$	${670}_{-31}^{+31}$	$141.18/172$	Obs4	16	${0.5585}_{-0.0063}^{+0.0065}$	${778}_{-12}^{+23}$	$145.32/169$
Obs3	17	${0.5322}_{-0.0052}^{+0.0053}$	${697}_{-31}^{+31}$	$202.26/172$	Obs4	17	${0.5606}_{-0.0060}^{+0.0061}$	${792}_{-23}^{+24}$	$200.99/168$
Obs3	18	${0.5351}_{-0.0049}^{+0.0050}$	${719}_{-31}^{+31}$	$170.80/171$	Obs4	18	${0.5631}_{-0.0057}^{+0.0059}$	${800}_{-20}^{+12}$	$163.63/170$
Obs3	19	${0.5407}_{-0.0047}^{+0.0048}$	${723}_{-31}^{+30}$	$173.93/171$	Obs4	19	${0.5624}_{-0.0056}^{+0.0057}$	${831}_{-23}^{+23}$	$203.28/169$
Obs3	20	${0.5360}_{-0.0047}^{+0.0048}$	${752}_{-31}^{+31}$	$173.03/171$	Obs4	20	${0.5659}_{-0.0056}^{+0.0058}$	${814}_{-23}^{+21}$	$158.42/170$

Download table as:  ASCIITypeset images: 1 2

Table 3.  Spectral Parameters of DISKBB and NTHCOMP

Component	Parameter	Obs1	Obs2	Obs3	Obs4
DISKBB	$k{T}_{0}$ (keV)	${0.277}_{-0.007}^{+0.007}$	${0.275}_{-0.005}^{+0.005}$	${0.233}_{-0.007}^{+0.007}$	${0.287}_{-0.006}^{+0.006}$
	Nor	${9441}_{-738}^{+842}$	${10125}_{-526}^{+576}$	${16454}_{-1550}^{+1765}$	${8126}_{-514}^{+551}$
NTHCOMP	$\mathrm{NTHCOMP}({\rm{\Gamma }})$	${1.483}_{-0.057}^{+0.054}$	${1.549}_{-0.026}^{+0.025}$	${1.727}_{-0.024}^{+0.021}$	${1.484}_{-0.039}^{+0.036}$
	${kT}_{e}\,({\rm{k}}{\rm{e}}{\rm{V}})$	${2.196}_{-0.071}^{+0.076}$	${2.334}_{-0.037}^{+0.039}$	${2.751}_{-0.071}^{+0.073}$	${2.215}_{-0.047}^{+0.049}$
	Nor	${0.113}_{-0.019}^{+0.020}$	${0.198}_{-0.016}^{+0.017}$	${0.318}_{-0.023}^{+0.023}$	${0.134}_{-0.016}^{+0.017}$
${\chi }^{2}/\mathrm{dof}$		3100.6/3144	3874.6/3638	3729.8/3417	3546.7/3387

Note.

^an_H is fixed at 0.36 × 10²² cm⁻².

Download table as:  ASCIITypeset image

Because the count rate is proportional to temperature to the power of 4 and area ($\mathrm{count}\ \mathrm{rate}\propto L\propto {{AT}}^{4}$). We investigated the correlations between the parameters of BBODYRAD and the QPO profile. Previous studies gave two contradictory results. Stiele et al. (2016) concluded that the variation in area dominates the flux oscillation of the mHz QPO, whereas Strohmayer et al. (2018) showed that variation is mainly owing to the temperature modulation. To study further which factor is more important to the mHz QPO, the Pearson correlation coefficient, r, was applied to investigate the linear relationship between two variables (i.e., count rate versus area and count rate versus temperature). The corresponding p-values, under the null hypothesis that two variables follow the Gaussian distribution with zero correlation coefficient, can be used to verify whether the spectral parameter has a strong linear relation to the oscillation count rates. A smaller p-value indicates a stronger correlation between two variables.

Table 4 lists the correlation coefficients and the corresponding p-values. For the Obs1 spectra, there is a strong correlation between T_bb and count rates with a linear correlation coefficient r = 0.84 and p-value of p = 3 × 10⁻⁶, but almost no correlation between the apparent area and count rates (r = 0.19, p = 0.41). This indicates that the oscillation is primarily attributed to the variation in temperature in an approximately constant apparent area (also refer to Figure 4). For the Obs2 spectra, the marginally positive correlation between T_bb and the count rates is r = 0.48 (p = 0.03), and the apparent area is highly correlated with the count rates with r = 0.91 (p = 2 × 10⁻⁸). For the Obs3 and Obs4 spectra, in addition to the marginally positive correlations between T_bb and the count rates, where r = 0.7 (p = 5 × 10⁻³) and r = 0.66 (p = 1 × 10⁻³) for Obs3 and Obs4, respectively, the apparent area is significantly correlated with the count rates, with r = 0.82 (p = 8 × 10⁻⁶) and r = 0.88 (p = 2 × 10⁻⁷), respectively. These results show that the variation in area dominates the mHz QPO with a constant temperature in Obs2, Obs3, and Obs4 but the variation in temperature dominates the mHz QPO with constant area in Obs1. Further implications of these results will be discussed in Section 4.

Table 4.  Correlation between Count Rate and Temperature (T) and between Count Rate and Area (A)

Note	T: r-value	T: p-value	A:r-value	A:p-value
Obs1	0.84	3 × 10−6	0.19	0.41
Obs2	0.48	0.03	0.91	2 × 10−8
Obs3	0.70	0.005	0.82	8 × 10−6
Obs4	0.66	0.001	0.88	2 × 10−7

Download table as:  ASCIITypeset image

We have utilized the HHT to characterize the HHT-based timing properties, extracted the instantaneous phase of the ∼8 mHz QPOs in 4U 1636-53, and constructed the modulation profiles and phase-resolved spectra for the complete cycle. The mHz QPO is a nonstationary periodic signal, for which the frequency or the amplitude would change with time. Because the dynamic power spectrum only yielded the mHz QPO frequency in a time interval, it cannot precisely derive the instantaneous frequency. Conversely, the HHT technique enables us to determine the instantaneous frequency for each data point, which provides an alternative point of view from which to study the nature of the mHz QPOs. More important is that the HHT may resolve the instantaneous phase, which is believed to be closely related to the physics behind the oscillations rather than the time. With the instantaneous phase, phase-resolved analyses can be performed (e.g., Hu et al. 2014).

In this work, we extracted the instantaneous phase of mHz QPO of 4U 1636-53. Variations in spectral parameters can be resolved for the whole mHz oscillation cycle. Our phase-resolved analyses that employ the HHT in this work have a phase resolution of 0.05 cycles for a complete QPO cycle. This shows that the HHT has the ability to provide more details in the variations of the spectral parameters. In a theoretical model, a mHz QPO is attributed to the variation of nuclear reaction rate modulations in the hot spot (Heger et al. 2007) that may cause the temperature oscillation or apparent variation in area. Recently, there have been two contradictory results for mHz QPO observations. Strohmayer et al. (2018) revealed significant oscillations at a frequency of ∼8 mHz in GS 1826-238 from observation with the Neutron Star Interior Composition Explorer (NICER) on 2017 September 9. This ∼8 mHz frequency and its amplitude are consistent with other accreting neutron-star systems with mHz QPO detections. Their phase-resolved spectra show that the mHz oscillation is produced by modulation of the temperature component of blackbody emission from the neutron-star surface with a constant emission area through the oscillation cycle. However, Stiele et al. (2016) provided a different conclusion for the mHz QPO of 4U 1636-53. Their phase-resolved spectra showed that mHz QPOs are a result of the periodically changing size of the hot spot with a constant temperature. In this work, the detailed and more precise phase-resolved spectra can exhibit the correlations shown in Section 3.2. The result reveals that the variation in area dominates the mHz QPO in Obs2, Obs3, and Obs4, but the variation in temperature has much stronger correlation with the oscillation count rate than the apparent area in Obs1. This indicates that variations in either temperature or apparent area of the hot spot may dominate the mHz QPO oscillation. The cause of the difference is not clear but it might be related to the spectral state of the source. Figure 1 in Lyu et al. (2015) showed the spectral states of these four observations in the color–color diagram. We find that the hard color of Obs1 is higher than the others. This is a possible clue to why Obs1 behaves differently than the other three. However, because there are only four observations with mHz QPO detections, more observations are required to verify whether the different behavior of Obs1 is due to the spectral state or just a coincidence.

One of the important factors in marginally stable burning is the accretion rate. The theoretical interpretation showed that the accretion rate is close to the Eddington accretion rate ($\sim {\dot{M}}_{\mathrm{Edd}}$) (Heger et al. 2007), but current observations indicated that the accretion rate is much lower ($\sim 0.1{\dot{M}}_{\mathrm{Edd}}$) for the sources with mHz QPOs. Fujimoto et al. (1981) suggested that the burst ignition of a neutron star depends on the accretion rate per unit area, $\dot{m}$, instead of the overall accretion rate. The local accretion rate ($\dot{m}$) does not need to be the same everywhere on the neutron-star surface (Bildsten 1998). The nuclear burning depends on the local accretion rate (Bildsten 1998; Heger et al. 2007). For a fast-spinning neutron star, a latitude change in the local accretion rate is influenced by the variation in the effective surface gravity from the equator to the pole (Strohmayer et al. 2018). Lyu et al. (2016) discovered that mHz QPOs are only associated with bursts with positive convexity. Because Maurer & Watts (2008) found that type-I X-ray bursts that ignite at the equator always have positive convexity, the mHz QPO is located on the equator of the neutron star. Strohmayer et al. (2018) roughly estimate the effective surface gravity of the neutron star in 4U 1636-53; the value on the pole can be ∼11% stronger that on the equator. They found that this value is sufficient to influence marginally stable burning by affecting the local accretion rate, and marginally stable burning on the equatorial belt on the neutron-star surface is possible (Altamirano et al. 2008; Strohmayer et al. 2018). Even though using the local accretion rate can explain marginally stable nuclear burning, a more complete theoretical model is required to explain the phenomena that we observed. Heger et al. (2007) calculated that marginally stable burning is one-dimensional; thus, the burning area is not explored (Strohmayer et al. 2018).

Data collected by XMM-Newton were analyzed in this work. The mHz QPOs are typically detected in low-energy bands (1–5 keV); thus, the energy range of XMM-Newton (0.1 to 15 keV) is beneficial for observing the mHz QPOs. XMM-Newton is a high-altitude satellite with an orbital period of 2789.6 minutes, which enables it to observe a source for a long time without interruption. Because the HHT prefers continuous data, the data collected by XMM-Newton are beneficial for using the HHT to analyze mHz QPOs. The high spectral resolution of XMM-Newton is also advantageous to the spectral analysis. In addition to XMM-Newton, the Nuclear Spectroscopic Telescope Array mission (NuSTAR) and NICER are possible choices that can be used to study mHz QPOs. Unfortunately, because NuSTAR focuses on high-energy X-ray bands (3–79 keV), it may not be suitable for studying mHz QPOs (typically detected in lower energy bands). The time resolution of NICER is 100 ns, and its spectral band overlaps those of RXTE and XMM-Newton (0.2–12 keV). The high time resolution enables us to analyze more timing properties by the HHT. NICER has a large effective area, which is beneficial for collecting photons. More photons have more timing information; thus, NICER would be a reasonable choice to perform timing analysis. However, NICER is a facility on board the International Space Station (ISS), and the observations are usually interrupted by Earth occultation for each ISS orbital cycle (∼90 minutes). On the other hand, X-ray telescopes that are expected to operate in the next decade, such as the enhanced X-ray Timing and Polarimetry mission (eXTP) and X-ray Imaging and Spectroscopy Mission (XRISM), although they have better time and spectral resolution in the low-energy band (∼0.5–10 keV), will be low-altitude satellites and the Earth occultation would induce too much alias in the HHT analysis, similar to NICER. These factors make the data collected by XMM-Newton the best data for analysis by the HHT for investigations of mHz QPOs. In the future, more XMM-Newton observations will enable us to investigate the physical causes of mHz QPOs by using the HHT to analyze more timing properties and spectral modulation of mHz QPOs.

We would like to thank an anonymous reviewer for helpful and valuable comments on the manuscript. We would like to thank Dr. Yi-Hao Su for useful advice regarding the HHT analysis. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). This work has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA's Goddard Space Flight Center. The HHT codes were provided by the Research Center for Adaptive Data Analysis in the National Central University of Taiwan. This research was supported by the grants MOST 108-2112-M-008-005- and 109-2112-M-008-004- from the Ministry of Science and Technology of Taiwan.

Facilities: ADS - , HEASARC - , XMM-Newton. -

Software: Science Analysis System (SAS) v.16.0.0. (Gabriel 2017), XSPEC (Arnaud 1996), astropy (Astropy Collaboration et al. 2013).