An estimation of the white dwarf mass in the Dwarf Nova GK Persei with NuSTAR observations of two states

Abstract

1 INTRODUCTION

Cataclysmic variables (CVs) are close binary systems consisting of a mass-accreting white-dwarf (WD) primary and a mass-donating companion. Gas overflowing from the Roche lobe of the companion accretes on to the WD surface, where gravitational energy of the gas is converted mainly into X-ray emission. CVs hosting a magnetized WD are further classified into ‘Polars’ and ‘Intermediate Polars’ (IPs), in which the WDs have magnetic field strengths of B ∼ 10^7 − 9 G and B ∼ 10^5 − 7 G, respectively.

An X-ray spectrum from an IP is a particular superposition of optically thin thermal emissions of various temperatures, from T_s downwards (Cropper, Ramsay & Wu 1998). To determine T_s, it is hence important to accurately measure both the hard X-ray continuum (e.g. Suleimanov, Revnivtsev & Ritter 2005; Yuasa et al. 2010), and the ratio of Fe XXV and XXVI lines at ∼7 keV (Fujimoto & Ishida 1997). This is because the former is sensitive to the hottest components (with temperature ∼T_s), whereas the latter tells us contributions from cooler components arising closer to the WD surface.

GK Persei, at an estimated distance of |$477 ^{+28}_{-25}$| pc (Harrison et al. 2013a), interestingly exhibits three distinct aspects of CVs; it behaves as an IP, as a Dwarf Nova, and exhibited a classical Nova explosion in 1901 (Hale 1901; Williams 1901). It repeats Dwarf Nova outbursts every 2–3 yr, each lasting for 2 months (e.g. Šimon 2002). During outbursts, the optical and X-ray luminosities both increase by a factor of 10–20.

A recent outburst from GK Persei started in 2015 March, and continued for 2 months (Wilber et al. 2015). During this outburst, Zemko et al. (2017) triggered a Target of Opportunity (ToO) observation with NuSTAR and measured a high-spin modulation even in a hard X-ray range. Suleimanov et al. (2016) also analysed the ToO data and constrained the WD mass as M_WD = 0.86 ± 0.02 M_⊙. The onset of this outburst was serendipitously caught by Suzaku, and the obtained data allowed Yuasa, Hayashi & Ishida (2016) to study the accretion geometry at the beginning of the outburst.

With NuSTAR, we observed GK Persei again, after the object returned to its quiescence. Although previous observations of the object in quiescence were unable to detect the hard X-ray component, the high sensitivity of NuSTAR has for the first time allowed us to detect its hard X-rays (typically in energies above ∼20 keV) in quiescence. This paper describes a combined analysis of the outburst and quiescence data from NuSTAR, and presents a new method to determine R_in, M_WD, and B of the WD in GK Persei utilizing the large change in |$\dot{M}$|⁠.

2 OBSERVATIONS AND DATA REDUCTION

The 2015 outburst of GK Persei started on 2015 March 6.84 ut (Wilber et al. 2015). As shown in Fig. 1, the ToO observation with NuSTAR (Harrison et al. 2013b) was conducted in the middle of the outburst from 2015 April 4 02:46:07 to April 6 15:10:35. The net exposures of focal plane module A and B (FPMA and B) are 42 ks each. The second NuSTAR observation, in quiescence, was performed from 2015 September 8 15:46:08 to September 11 02:04:09 with a net exposures of 72 ks. The log of the two observations is given in Table 1.

Figure 1.

Optical light curves of GK Persei without a filter from AAVSO International Database, and the Swift/BAT (Krimm et al. 2013) 15–150 keV count rate history from Swift/BAT X-ray Transient Monitor web site. The shaded regions indicate the two NuSTAR observations.

Open in new tabDownload slide

We utilized the data analysis software package HEASOFT version 6.20 and a detector calibration data base NuSTAR CALDB version 20170222, both released and maintained by HEASARC at NASA Goddard Space Flight Center. Photon events in the data sets were extracted with an exclusive data reduction software for NuSTAR ‘nupipleine’ version 0.4.6 and ‘nuproducts’ version 0.3.0. The on-source events were accumulated from a circular region with a radius of 150 arcsec (in outburst) and 80 arcsec (in quiescence) centred on the source. The background data were accumulated over a region outside a circle of radius of 170 and 100 arcsec in outburst and quiescence, respectively. The X-ray spectra were analysed and fitted with XSPEC version 12.9.1 (Arnaud 1996).

Generally, the X-ray emission from an IP is pulsed at its P. In fact, pulsations of GK Persei in the X-ray band have been detected at P = 351 s both in outbursts and quiescence (e.g. Watson et al. 1985; Norton, Watson & King 1988; Ishida et al. 1992). In the 2015 outburst observation, Zemko et al. (2017) clearly detected the 351 s pulsation both in the 3–10 and 10–79 keV ranges. In the quiescence observation by NuSTAR, we detected a faint pulsation with a modulation amplitude of ∼10 per cent in the 3–50 keV range. In this paper, we concentrate on spectral analysis and postpone the study of this pulsation for the next publication.

3 ANALYSIS AND RESULTS

Fig. 2 shows 3–50 keV spectra of the outburst and quiescence observations. The background has been subtracted, but the instrumental response has not been removed. Data of FPMA and FPMB are separately plotted. As reported by Zemko et al. (2017) and Suleimanov et al. (2016), the hard X-ray continuum is detected up to 70 keV during the outburst. In quiescence, the source is detected up to 50 keV for the first time. The 3–50 keV count rate of FPFA plus FPMB was 18.09 ± 0.02 count s⁻¹ in outburst, and 1.080 ± 0.006 count s⁻¹ in quiescence. Thus, the outburst data have 17.5 times higher count rate than those in quiescence. The spectra, particularly the outburst data, exhibit Fe–K line complex at ∼6.4 keV. From this energy, we regard the lines as mainly of fluorescence origin (from the WD surface and/or the accreting cold matter), rather than ionized lines from the accretion columns.

Figure 2.

(a) Time-averaged FPMA (black) and FPMB (red) spectra of the outburst (brighter) and quiescence (dimmer) observations. (b) Ratios of the FPMA plus FPMB spectra between the two observations.

Open in new tabDownload slide

In the bottom panel, we show the ratio between the two spectra. It reveals three features of the outburst spectrum, in comparison with that in quiescence. Namely, a stronger low-energy absorption, the stronger Fe–K line, and a continuum break at ∼20 keV.

With the model thus constructed, we first fitted the outburst spectrum in the 5–50 keV range, because the cemekl model is not available above 50 keV. A partial covering absorption model was applied to the spectral model in addition to a single column absorber. As shown in Fig. 3(a), this model approximately reproduced the spectrum, but the fit was formally not acceptable under a 90 per cent confidence level, with the reduced chi-squared of χ²/ν = 2.26 for 133 degrees of freedom even including 1 per cent systematic error. In fact, as shown in Fig. 3(b), significant residuals were seen in the low-energy band (<10 keV).

Figure 3.

(a) The NuSTAR FPMA (black) and FPMB (red) spectra in outburst, and their best-fitting models with a partial covering absorption. The thermal component, reflection component and Fe–K line are indicated by the dash–dotted, dash–dot-dotted, and dashed lines, respectively. (b) The residuals from the 5–50 keV fit in panel (a). (c) The residuals of the 15–50 keV fit. (d) The quiescence spectrum and its best-fitting model. (e) The residuals of the fit in panel (d).

Open in new tabDownload slide

The above fit failure to the outburst spectra is not surprising, because X-ray spectra of IPs are often subject to strong and complex absorption that is not modelled by partial absorption (e.g. Ezuka & Ishida 1999). Since refining the absorption model is beyond the scope of this paper due to lack of constraining data other than the continuum shape, we have resorted to discarding low-energy ranges until the effects of complex absorption become negligible (see Ezuka & Ishida 1999, for a similar method utilized to avoid complex absorption from affecting the spectral fitting result). By limiting the fit range to 15–50 keV, the fit to the outburst spectrum has become acceptable even with a single column density absorption. The range of T_s constrained in this way, T_s = 19.4 ± 0.8 keV, approximately accommodates the value of T_s = 17.9 keV obtained using the 5–50 keV range (though the fit was unacceptable).

The 5–50 keV quiescence spectrum has been reproduced successfully by the spectral model with a single column absorber. Therefore, we have finally conducted a simultaneous fitting using the 15–50 keV band of the outburst spectrum and the 5–50 keV band of the quiescence spectrum. The inclination angle and the abundance were set in common, while the other parameters were allowed to vary independently. This fitting including 1 per cent systematic error has become acceptable with χ²/ν = 1.01 for 191 degrees of freedom. The fit result and the best-fitting parameters are presented in Fig. 3 and Table 2, respectively. Errors are at 90 per cent confidence level. The shock temperature was constrained as |$T_{\rm s} = 19.7^{+1.3}_{-1.0}$| keV in outburst, and |$T_{\rm s} = 36.2^{+3.5}_{-3.2}$| keV in quiescence; thus, T_s was significantly higher in the latter. The 15–50 keV absorbed flux in outburst was 33 times higher than that in quiescence. The 15–50 keV luminosity is also derived as |$1.2^{+0.2}_{-0.3} \times 10^{34}$| erg s⁻¹ in outburst, and |$3.5^{+0.3}_{-0.5} \times 10^{32}$| erg s⁻¹ in quiescence, with the estimated distance of 477 pc (Harrison et al. 2013a).

For our purpose, we need to calculate the total X-ray flux F which is thought to represent the gravitational energy release from R_in to R_WD. Starting from the absorbed 15–50 keV flux, F was derived in the following way. First, the absorption and the reflection were removed. Second, the flux above 50 keV was included by extrapolating the best-fitting model up to 100 keV. Finally, the contribution below 15 keV was incorporated by integrating the best-fitting model down to 0.01 keV. The flux above 100 keV and below 0.01 keV are both estimated to be ≪0.01F. The difference in F between the two spectra amounts to a factor of 65.

Since our final fit to the outburst spectrum was obtained by discarding the data below 15 keV, obviously no information was obtained on the Fe–K line. Hence, we fitted the 5–9 keV spectra of the outburst and quiescence to constrain the Fe–K line equivalent width (EW). With a simple model of a single-temperature bremsstrahlung, a Gaussian and a single column absorption model, the EW was measured to be 192 ± 13 eV and |$52^{+34}_{-26}$| eV in outburst and quiescence, respectively, as presented in Table 2. The latter is consistent with the value obtained by the 5–50 keV simultaneous fitting.

4 ESTIMATION OF THE WD MASS

Figure 4.

Contours (dotted orange curves) of T_s given by equation (6), shown on the M_WD/M_⊙–R_in/R_WD plane. Allowed regions for T_s in outburst and quiescence are indicated by blue and green, respectively. The dashed red line (vertical) shows a value of M_WD that satisfies γ = 3.9 between the two measurements.

Open in new tabDownload slide

We can now determine the value of M_WD in Fig. 4 so that the two R_in values satisfy equation (11) as indicated by a pair of horizontal lines in Fig. 4; |$R_{\rm in}^{\rm out} = 1.9^{+0.4}_{-0.2}\ R_{\rm WD}$|⁠, |$R_{\rm in}^{\rm qui} = 7.4^{+2.1}_{-1.2}\ R_{\rm WD}$| and M_WD = 0.87 ± 0.05 M_⊙ are yielded with R_WD = (6.6 ± 0.4) × 10⁸ cm. The obtained |$R_{\rm in}^{\rm qui}$| and |$R_{\rm in}^{\rm out}$| satisfy the accretion condition R_in < R_Ω ∼ 11 R_WD with P = 351 s, M_WD = 0.87 M_⊙ and equation (3). Thus, the clear increase in T_s, observed in the transition from outburst to quiescence (Table 2), has been successfully explained as a factor of γ ∼ 4 increase in R_in, in response to the decrease in |$\dot{M}$| by a factor of 120: |$\dot{M}$| is obtained from equation (9) and F as |$\dot{M}^{\rm out} = (1.2 \pm 0.4)\times 10^{18}\ {\rm g\ s^{-1}}$| and |$\dot{M}^{\rm qui}\ =\ (1.0 \ \pm \ 0.3) \times 10^{16}\ {\rm g\ s^{-1}}$|⁠, respectively. An advantage of this method is that any systematic uncertainty involved in the coefficient of equation (7) cancels by taking the ratio of the two equations (equation 10). Further discussion continues in Section 5.4.3.

5 DISCUSSION

5.1 Comparison between the two observations

In this paper, we analysed a pair of NuSTAR spectra of GK Persei acquired in 2015. In the 5 months from the outburst observation on April 4 to the quiescence one on September 8, the optical emission from GK Persei diminished by ∼3.2 magnitude or ∼19 times (Fig. 1). Meanwhile, the 3–50 keV FPMA+FPMB count rate decreased by a factor of 17, and the absorbed 15–50 keV flux by a factor of 33. The two factors become different because the outburst spectrum is more absorbed (Fig. 2), and hence the count rate which is more weighted towards lower energies changed less than that of the flux which is more weighted towards higher energies. Correcting these spectra for the respective absorption, and extrapolating the best-fitting models to >50 keV and <15 keV, the 0.01–100 keV unabsorbed total X-ray flux is inferred to have changed by 65 times.

In addition to the changes in the X-ray flux and absorption, we detected a clear increase in T_s from the outburst (∼20 keV) to the quiescence (∼36 keV) observations. Employing the disc–magnetosphere interaction model of Ghosh & Lamb (1979), the change in T_s has successfully been interpreted as due to a factor of ∼4 change in R_in, in a negative correlation with |$\dot{M}$|⁠. Considering that the gravitational potential drop available for the X-ray emission (from R_in to R_WD) thus became deeper in quiescence, the total X-ray luminosity change has been converted to a factor of 120 difference in |$\dot{M}$|⁠. For reference, the temperature change we observed is qualitatively consistent with the report by Zemko et al. (2017), that the Swift/XRT light curve of the hardness ratio indicated a temperature decrease as the outburst proceeded towards its peak, and the very high value of T_s measured by Yuasa et al. (2016) at the outburst onset.

As presented in Table 2, the EW of the Fe–K line (nearly neutral component) was ∼4 times higher in outburst. The line is usually ascribed to two emitting sources: the ambient matter as represented by |$\dot{M}$| and N_H, and the WD surface as represented by reflection. Evidently, both |$\dot{M}$| and N_H were higher in outburst, so that the Fe–K line EW from the first source must be higher as well. Furthermore, as discussed later in Section 5.4.2, the standing shock is considered to come closer to the surface when |$\dot{M}$| increases, because the higher density would increase the volume emissivity in the accretion columns and higher shock temperature can be dissipated within the boundary conditions. This will in turn increase the solid angle of reflection, and yield a high EW from the second source. The higher EW observed in outburst may be explained qualitatively as a combination of these two effects.

5.2 Comparison with previous optical results

5.3 Comparison with previous X-ray results

Let us revisit the past X-ray result with the PCA and HEXTE onboard RXTE, obtained during an outburst by Suleimanov et al. (2005). They measured T_s = 21 ± 3 keV, and derived M_WD = 0.59 ± 0.05 M_⊙ assuming R_in ≫ R_WD (the author noted that M_WD would be underestimated). At that time, the total X-ray flux in the 0.1–100.0 keV range was measured to be 8.86 × 10⁻¹⁰ erg cm⁻² s⁻¹. When we use the present cemekl model of the same T_s, the 0.01–100 keV flux is re-estimated as 1.0 × 10⁻⁹ erg cm⁻² s⁻¹, which falls in between the present two measurements (Table 2). Then, compared with them, the value of R_in during the RXTE observation is estimated, from equation (8), as |$R_{\rm in} \simeq 1.5\ R_{\rm in}^{\rm out} \simeq 0.37\ R_{\rm in}^{\rm qui} \simeq 2.8\ R_{\rm WD}$|⁠. Substituting this value and T_s = 21 keV into equation (6), or equivalently referring to Fig. 4, M_WD ∼ 0.8 M_⊙ is derived. This revised mass is probably consistent with our result when various errors are taken in account. Our result is also consistent with the mass estimation M_WD = 0.90 ± 0.12 M_⊙ by Brunschweiger et al. (2009) with Swift/BAT within errors, as already referred to in Section 1.

The present outburst data were already analysed by Zemko et al. (2017). They derived |$T_{\rm s} = 16.2^{+0.5}_{-0.4}\ {\rm keV}$| by the 3–50 keV broad-band fitting, wherein the NuSTAR data are combined with those from the Chandra MEG and HETG. Their T_s value is ∼18 per cent lower than our outburst result. This discrepancy may be caused by difference of emission models. The mkcflow model they employed is a superposition of the mekal thermal emission model, like the cemekl model we used, but the emission measure of mkcflow is weighted by the inverse of the bolometric luminosity at each temperature T. Because the bolometric flux is ∝T^1/2 when only the bremsstrahlung continuum is considered, the differential emission measure becomes ∝(T/T_s)^−1/2, and α = 0.5 by equation (4). It is slightly different from that of cemekl, α = 0.43, and will make the composite spectrum more weighted towards higher temperatures. In the fit by Zemko et al. (2017), this effect is considered to be compensated by the lower T_s. (In the relevant temperature range, α would not change very much even considering the lines.)

Suleimanov et al. (2016) also analysed the same outburst data and obtained M_WD = 0.86 ± 0.02 M_WD. This is fully consistent with our estimate. In deriving this result, however, they employed a method that differs from ours in two points: T_s and R_in. They fitted the 20–70 keV spectrum in outburst with their newly calculated spectral model (PSR model), and obtained T_s ∼ 26.3 keV. They also employed, for an illustrative purpose, a single temperature bremsstrahlung model and obtained its temperature as 16.7 ± 0.2 keV; via equation (2) in Suleimanov et al. (2016), this was converted to a consistent shock temperature of T_s = 26.0 ± 0.3 keV. For consistency, we fitted the 20–70 keV spectrum with the bremsstrahlung model, and obtained the temperature as 16.6 ± 0.3 keV. Therefore, the present data analysis is consistent with theirs. They also derived R_in = 2.8 ± 0.2 R_WD by the power density spectral analysis (Revnivtsev et al. 2009). After all, their T_s is 1.3 times higher than our T_s, and their R_in is 1.5 times larger than ours. These differences in T_s and R_in happened to cancel out, to yield the two M_WD estimates which are very close to each other.

5.4 Systematic uncertainties of the mass estimation

So far, we considered only statistic errors. Here, let us evaluate possible systematic errors that can affect our result.

5.4.1 Emission models of the accretion column

In this paper, we assumed the accretion columns to have a cylindrical shape, and hence employed the cemekl model with α = 0.43 in equation (4). Recently, emission models with a dipole geometry for the accretion column have been developed, including ACRAD model by Hayashi & Ishida (2014a,b) and the PSR model by Suleimanov et al. (2016). We thus refitted the spectra with the PSR model (ipolar model in XSPEC). The energy range below 7 keV in quiescence was ignored in this analysis because ipolar model has abundances fixed to 1 Z_⊙ and cannot reproduce the Fe emission lines which are well described with sub-solar Fe abundance. The shock temperature was then obtained as 20.4 ± 0.6 keV in outburst and |$37.3 ^{+3.9}_{-3.3}$| keV in quiescence, and the WD mass was constrained as M_WD = 0.88 ± 0.05 M_⊙. All these values agree well with the results in Section 3 within the statistic errors. Therefore, we consider that slight differences of the emission models, namely the detailed morphological and emissivity structures of the post-shock region, have insignificant impact on the mass estimation method presented above, at least, when applied to GK Persei.

5.4.2 Shock height

As presented in equation (5), T_s depends on the shock height h, which is thought to negatively correlate with |$\dot{M}$|⁠. This effect was theoretically calculated and incorporated in equation (6). However, even if the effect of h is ignored (i.e. h = 0), the value of M_WD changes by less than 0.5 per cent, which is much smaller than the statistic error. Therefore, the value of h itself would not affect the WD mass estimation in GK Persei.

In our spectral analysis, the reflection model was included to represent the reflection effect on the WD surface. The solid angle of the reflection was fixed to 2π simply assuming that the shock heating occurs just above the WD surface. In reality, h are calculated as 0.014 R_WD and 0.026 R_WD by equation (5), corresponding to the solid angle of 1.8π and 1.7π in outburst and quiescence, respectively. We thus repeated the spectral fitting using these solid angles, to find that neither |$T_{\rm s}^{\rm out}$| nor |$T_{\rm s}^{\rm qui}$| changes by more than 0.2 per cent Therefore, the result of M_WD is not affected either.

5.4.3 Alfven radius

As expressed by equation (7), we assumed that R_in is equal to the Alfvén radius. This formalism by Elsner & Lamb (1977), which considers spherical accretion, includes two possible uncertainties. One is the coefficient of equation (7). In the case of disc accretion in rotating magnetic neutron star systems, Ghosh & Lamb (1979) argued that the actual inner-disc radius is almost half the Alfvén radius. Therefore, the estimated values of R_in are probably subject to an uncertainty by a constant factor. However, this uncertainty cancels out in our work by taking the ratio |$\gamma = R_{\rm in}^{\rm qui}/R_{\rm in}^{\rm out}$|⁠. Hence, our M_WD estimation is free from the uncertainty, because in Fig. 4 we utilize γ rather than the actual values of R_in.

The other uncertainty is the power-law index Γ of |$\dot{M}$| in equations (7), (8) and (10), for which we employ Γ = −2/7. A recent 3D simulation of accretion flows around a neutron star yielded Γ = −1/5 (Kulkarni & Romanova 2013). In addition, Suleimanov et al. (2016) observationally evaluated |$\Gamma = -0.2^{+1.0}_{-1.5}$|⁠. We thus quote a systematic uncertainty by ΔΓ ∼ 0.1, which translates into ∼7 per cent systematic errors in M_WD.

Adding up all these uncertainties, the overall systematic error in the WD mass estimate becomes comparable to the statistical error. Including this in quadrature, we quote our final mass determination as M_WD = 0.87 ± 0.08 M_⊙, with R_WD = (6.6 ± 0.6) × 10⁸ cm.

5.5 Strength of magnetic field on WD surface

Since M_WD, R_WD and R_in were determined, the strength of the magnetic field on the WD surface can be now estimated to be B ∼ 5 × 10⁵ G using equation (7) and the distance of the object 477 pc (Harrison et al. 2013a). This value is consistent with the typical magnetic field strength of IPs. However, unlike M_WD, this result is subject to the uncertainties discussed in Section 5.4.3. If the coefficient of equation (7) has an uncertainty by a factor of 2 for example, B changes by a factor of 3–4.

Magnetic field measurements of Polars and strongly magnetized IPs (∼10⁷ G) have been made by detecting spin-modulated polarization in the near-ultraviolet to near-infrared bands (e.g. Piirola et al. 2008). However, the magnetic field of ∼10⁵ G on the WD surface cannot be measured at present. Therefore, the present method provides the only way to measure relatively weak magnetic field of IPs, even though it has a relatively poor accuracy due to the above-mentioned model uncertainty.

6 CONCLUSION

Analysing the outburst and quiescence data of GK Per obtained with NuSTAR, we found that the 0.01–100 keV unabsorbed flux was |$3.6^{+0.5}_{-0.8} \times 10^{-9}\ {\rm erg}\ {\rm s}^{-1}\ {\rm cm}^{-2}$| and 5.5|$^{+0.5}_{-0.9} \times 10^{-11}\ {\rm erg}\ {\rm s}^{-1}\ {\rm cm}^{-2}$|⁠, respectively, with a factor of 65 difference. Analysing the 5–50 keV or 15–50 keV spectra using a multitemperature spectral model, the shock temperature was determined as |$19.7^{+1.3}_{-1.0}$| keV in outburst, and |$36.2^{+3.5}_{-3.2}$| keV in quiescence. Assuming that this temperature difference is caused by a compression of the magnetosphere and the associated change in R_in, we determined the WD mass in GK Per as M_WD = 0.87 ± 0.08 M_⊙, together with the radius as R_WD = (6.6 ± 0.6) × 10⁸ cm (including a 7 per cent systematic error). The values of R_in, relative to R_WD, were derived as |$R_{\rm in}/R_{\rm WD}=1.9^{+0.4}_{-0.2}$| in outburst and |$R_{\rm in}/R_{\rm WD}=7.4^{+2.1}_{-1.2}$| in quiescence, and the mass accretion rate is estimated to have changed by a factor of 120 between the two observations. Combined with optical observations, the inclination angle of GK Per is tightly constrained as 63° ≤ i ≤ 73°. We also estimated the magnetic field of the WD as ∼5 × 10⁵ G although it is subject to large systematic error uncertainty. The overall results demonstrate the power of our mass determination method using X-ray luminosity changes, wherein some major systematic uncertainties cancel out.

Acknowledgements

This work has made use of data obtained with the NuSTAR mission, and software obtained from the High Energy Astrophysics Science Archive Research Center at NASA Goddard Space Center. The authors also acknowledge the use of X-ray monitoring data provided by Swift data archive, and optical data from AAVSO International Database contributed by observers worldwide. YW is supported by the Junior Research Associate Programme in RIKEN.

REFERENCES