NuSTAR and XMM-Newton observations of SXP 59 during its 2017 giant outburst

TL;DR

This study analyzes the 2017 giant outburst of Be X-ray pulsar SXP 59 using XMM-Newton and NuSTAR data, revealing energy-dependent pulse profiles, a transition from super-critical to sub-critical accretion regimes, and complex spectral components.

Contribution

First detailed multi-epoch X-ray analysis of SXP 59 during its giant outburst, highlighting spectral and timing properties and implications for neutron star magnetic field and accretion physics.

Findings

Pulse fraction increases with energy, saturating above 10 keV.

Transition from double-peaked to single-peaked pulse profiles indicates regime change.

Presence of a hot blackbody component not explained by current accretion models.

Abstract

The Be X-ray pulsar (BeXRP) SXP 59 underwent a giant outburst in 2017 with a peak X-ray luminosity of $1.1\times10^{38}$ erg~s$^{-1}$. We report on the X-ray behaviour of SXP 59 with the XMM--Newton and NuSTAR observations collected at the outburst peak, decay, and the low luminosity states. The pulse profiles are energy dependent, the pulse fraction increases with the photon energy and saturates at $\sim$ 65% above 10 keV. It is difficult to constrain the change in the geometry of emitting region with the limited data. Nevertheless, because the pulse shape generally has a double-peaked profile at high luminosity and a single peak profile at low luminosity, we prefer the scenario that the source transited from the super-critical state to the sub-critical regime. This result would further imply that the neutron star (NS) in SXP 59 has a typical magnetic field. We confirm that the soft excess revealed below 2 keV is dominated by a cool thermal component. On the other hand, the NuSTAR spectra can be described as a combination of the non-thermal component from the accretion column, a hot blackbody emission, and an iron emission line. The temperature of the hot thermal component decreases with time, while its size remains constant ($R \sim 0.6$ km). The existence of the hot blackbody at high luminosity cannot be explained with the present accretion theories for BeXRPs. It means that either more sophisticated spectral models are required to describe the X-ray spectra of luminous BeXRPs, or there is non-dipole magnetic field close to the NS surface.

Tables

Table 1. Table 1: Log of XMM–Newton and NuSTAR observations.

Obs Date	Observatory	ObsID	Exposure	Epoch	$\nu$	$\dot{\nu }$
			(ksec)	(MJD-57850)	(${10}^{-2}$ Hz)	(${10}^{-11}$ Hz s-1)
2017 Apr 14	XMM–Newton	0740071301	14	7.617424	1.69633(4)	…
2017 Apr 12–13	NuSTAR	30361001002	70	5.786459	1.69570(2)	3.4(4)
2017 Apr 24–26	NuSTAR	50311001002	153	5.786930	1.69669(7)	2.02(7)
2017 Aug 12–13	NuSTAR	50311001004	82	127.093712	1.700745(7)	…

Table 2. Table 2: Spectra are fitted with the model of tbabs*(apec+bbodyrad+bbodyrad+cutoffpl+gaussian) . a Unabsorbed luminosity is calculated in 0.5–79 keV by assuming a distance of 62.1 kpc. b Unabsorbed luminosity is calculated in 3–79 keV. All errors are in 90% confidence level.

Parameters	XMM+NuSTAR	NuSTAR	NuSTAR	NuSTAR
	Apr 12–14	Apr 12–13	Apr 24–26	Aug 12–13
$nH$ (${10}^{22}$ cm-2)	${0.10}_{-0.03}^{+0.03}$	0.10 (fixed)	0.10 (fixed)	0.10 (fixed)
$k{T}_{\mathrm{apec}}$ (keV)	${1.00}_{-0.12}^{+0.16}$	–	–	–
Normapec ($\times {10}^{-4}$)	${5.3}_{-2.4}^{+2.2}$	–	–	…
$k{T}_{\mathrm{BB}}^{\mathrm{low}}$ (keV)	${0.19}_{-0.01}^{+0.01}$	–	–	…
Norm${}_{}{}^{\mathrm{low}}{}_{\mathrm{BB}}{}^{}$	${226.5}_{-109.2}^{+184.8}$	–	–	…
$k{T}_{\mathrm{BB}}^{\mathrm{high}}$ (keV)	${3.36}_{-0.25}^{+0.21}$	${4.09}_{-0.16}^{+0.14}$	${3.68}_{-0.14}^{+0.14}$	${1.66}_{-0.34}^{+0.52}$
Norm${}_{}{}^{\mathrm{high}}{}_{\mathrm{BB}}{}^{}$ ($\times {10}^{-3}$)	${9.5}_{-2.2}^{+2.2}$	${10.1}_{-1.9}^{+1.8}$	${9.7}_{-1.3}^{+1.3}$	${10.1}_{-7.1}^{+21.6}$
$\Gamma$	${0.64}_{-0.04}^{+0.05}$	${0.85}_{-0.07}^{+0.07}$	${1.11}_{-0.07}^{+0.08}$	${0.78}_{-0.49}^{+0.35}$
$E_{\mathrm{cut}}$ (keV)	${17.8}_{-1.0}^{+1.2}$	${22.7}_{-2.1}^{+2.5}$	${31.0}_{-4.2}^{+5.8}$	${11.7}_{-3.3}^{+6.3}$
Normcut-off ($\times {10}^{-3}$)	${2.8}_{-0.1}^{+0.1}$	${4.2}_{-0.3}^{+0.3}$	${3.7}_{-0.3}^{+0.3}$	${0.25}_{-0.14}^{+0.13}$
$E_{\mathrm{Gau}}$ (keV)	${6.30}_{-0.09}^{+0.09}$	${6.31}_{-0.10}^{+0.10}$	${6.26}_{-0.15}^{+0.19}$	–
$\sigma$ (keV)	${0.33}_{-0.09}^{+0.10}$	${0.31}_{-0.10}^{+0.11}$	${0.55}_{-0.19}^{+0.34}$	–
Normgauss ($\times {10}^{-5}$)	${3.9}_{-0.9}^{+1.0}$	${4.6}_{-1.2}^{+1.3}$	${3.8}_{-1.2}^{+1.9}$	–
$L_X$ (${10}^{38}$ erg s-1)	1.11a	1.05b	0.59b	0.032b
${\chi }^2/dof$	1806.2/1476	1047.8/1040	1119.3/975	408.0/403

Full text

1 Introduction

In general, Be X-ray binary (BeXRB) consists of a young neutron star (NS) orbiting a Be type star. Most BeXRBs are transient in X-ray, and their variabilities are often classified into two types of outbursts (see Bildsten et al., 1997; Reig, 2011, for reviews). Type I outbursts are less energetic ($L_{\rm peak}<10^{37}$ erg s*-1*) and occur regularly as the enhancement of accretion during the periastron passage. On the other hand, type II outbursts are rare and not fixed to the orbital phase. Their X-ray luminosity can exceed the Eddington luminosity for a NS. In particular, the peak X-ray luminosities of the 2016-17 outburst of SMC X-3 (e.g. Weng et al., 2017; Zhao et al., 2018) and 2017-18 outburst of Swift J0243.6+6124 (Doroshenko et al., 2018; Wilson-Hodge et al., 2018; Tao et al., 2019) are beyond $10^{39}$ erg s*-1*, that is the threshold of ultraluminous X-ray sources (ULXs, Kaaret et al., 2017).

Hundreds of high-mass X-ray binaries (HXMBs) have been detected in the Galaxy and the Magellanic Clouds, and more than half of them are BeXRBs (Liu et al., 2005, 2006). In addition to the short distance (62.1 kpc; Hilditch et al., 2005; Graczyk et al., 2014; Scowcroft et al., 2016), the Small Magellanic Cloud (SMC) has a large star formation rate (150 times of the Galaxy; Harris & Zaritsky, 2004) and low interstellar absorption (Zaritsky et al., 2002; Willingale et al., 2013). It thus provides an ideal and large sample of BeXBRs for a detailed study in multibands (e.g. Rajoelimanana et al., 2011; Bird et al., 2012; Coe & Kirk, 2015). Historically, many X-ray observatories (e.g. ROSAT, RXTE, Chandra, XMM–Newton) had spent a lot of time to survey HXMBs in the SMC (e.g. Kahabka et al., 1999; Galache et al., 2008; Sturm et al., 2013; Haberl & Sturm, 2016; Yang et al., 2017). Since 2016 June, the Neil Gehrels Swift Observatory has started a high cadence shallow (with a typical exposure of 60 s) survey of the SMC in order to monitor X-ray variabilities of BeXRBs by taking its advantage of rapid slewing. During the first year operation, the Swift SMC survey (S-CUBED) successfully detected the type II outbursts from SMC X-3, SXP 59, and SXP 6.85 (see Kennea et al., 2018, for more details).

SXP 59 was identified as an X-ray pulsar ($P=59.0\pm 0.2$ s) in 1998 due to its outburst, while the pulsations with the same period were also revealed in the ROSAT archive data (Marshall et al., 1998). The orbital period of $\sim 122.1$ d was reported with both the RXTE observations (Galache et al., 2008) and the OGLE I-band light curves (Bird et al., 2012). S-CUBED detected the onset of giant outburst from SXP 59 on 2017 March 30 (Kennea et al., 2017), and the Swift TOO observations were triggered to follow the outburst. The source reached a peak luminosity on 2017 April 7 ($\sim 4.6\times 10^{37}$ erg s*-1* in 0.5–10 keV), then exponentially declined with an time-scale of $\sim 15.9$ d, and returned to the pre-outburst flux level on 2017 June 6 (Kennea et al., 2018). Investigating the XMM–Newton TOO observation performed around the peak of outburst, La Palombara et al. (2018) revealed a soft excess below 2 keV in addition to the primary power-law component. Since the double-peaked pulse profile detected at the high-luminosity level, they also speculated that the source was at the super-critical state having a fan-beam emission geometry.

In this paper, we carry out a detailed analysis on the high-quality data obtained from the XMM–Newton and another three NuSTAR observations executed at different flux levels to explore the spectral evolution of SXP 59 during its 2017 giant outburst. Section 2 describes the observations together with the data analysis, and summarizes our results. We discuss the physical implications of these results in Section 3.

2 Data Analysis

2.1 X-ray observations

In 2017, the Neil Gehrels Swift Observatory carried out 92 observations on SXP 59, including 65 S-CUBED observations. In order to avoid the pile-up effect, the TOO observations around the outburst peak were executed with the window timing (WT) mode instead of the photon counting (PC) mode. The Swift/XRT data are processed with the packages and tools available in heasoft 6.24. The software xrtpipeline is used with standard quality cuts for the initial event cleaning. We extract the source light curves in 0.3–10 keV from a circle of 15 pixels centred at the source position, and the background light curves from an annulus region with the radii of 15 and 30 pixels. The source light curves are corrected for the telescope vignetting and point spread function losses with the task xrtlccorr, and then are subtracted by the scaled background count rate to generate the net light curves (Figure 1). When the source was in quiescence state, it can hardly be detected by the S-CUBED observations due to their short exposures. Following the work in Kennea et al. (2018), we adopt five counts as the threshold of detection, and calculate the upper limits for non-detections. Based on the following XMM–Newton spectral fitting results, we convert the count rate to the flux, and hence the luminosity assuming a distance of 62.1 kpc. The derived count rate to luminosity ratios for the WT and the PC modes are of 1 count s*-1* $\sim 3.55\times 10^{37}$ erg s*-1* and $\sim 3.77\times 10^{37}$ erg s*-1*, respectively. As can be seen in Figure 1, the giant outburst lasted for about two months with a fast-rise near exponential decay profile. These results are consistent with those reported in Kennea et al. (2018) (Figure 13 in their paper).

In this work, we analyze the XMM–Newton observation performed on 2017 April 14, which was free of the background contamination. The data collected with the XMM–Newton EPIC instrument are reduced using the Science Analysis System software (sas) version 14.0.0. Both the pn and MOS data were taken in small window mode in order to minimize the pile-up effect. We exclude all events at the edge of CCD and from bad pixels by setting FLAG=0, select the pn events with PATTERN in the 0-4 range, and the MOS data with PATTERN$\leq$12. The source photons are extracted from a circle aperture with a radius of 30 arcsec, and the background is taken from the same CCD chip as the source within a circle of radius 50 arcsec.

NuSTAR is the first direct-imaging hard X-ray telescope, consisting of two focusing instruments and two focal plane modules, i.e. Focal Plane Modules A and B (hereafter FPMA and FPMB; Harrison et al., 2013). There are three NuSTAR observations carried out at the outburst peak, decay, and the low luminosity states, respectively (Table 1 and Figure 1). The source events are extracted from circular region with the radius of 60–100 arcsec, depending on the count rates. Meanwhile, the background photons are extracted from the source-free region with a radius of 120″. These data are processed with the task nupipeline, the spectra and the light curves are produced with the command nuproducts. It is worth to note that, the first NuSTAR data were made 1–2 d before the XMM–Newton observation. That is, these two observations are quasi-simultaneous.

2.2 Spectral analysis

Both NuSTAR and XMM–Newton spectra are fitted by empirical models most often used in the literature with the heasoft X-ray spectral fitting package xspec (Arnaud, 1996). All models in this paper also include the interstellar absorption (tbabs in xspec). The NuSTAR spectra are grouped with grppha to ensure at least 30 counts per bin. The FPMA and the FPMB spectra are fitted simultaneously, with a constant multiplicative factor to compensate for calibration differences. The FPMA constant is fixed at unity, whilst that for the FPMB is allowed to vary, with the yielded values in the range of 1.00–1.06.

Because the NuSTAR and XMM–Newton data are operated in different energy ranges (3–79 keV and 0.5–10 keV), some emission component might be caught by only one of them. Thus, we firstly decompose the spectral components with NuSTAR and XMM–Newton spectra separately, and aim to achieve a common model for the broadband spectra. We begin by fitting the cut-off power-law component to the first NuSTAR observation. The derived reduced chi-square is of $1.09$ ($\chi^{2}/dof=1139.0/1045$), the iron line feature and residuals at low and high-energy bands are displayed in the top panel of Figure 2. Thus, a Gaussian line at $\sim 6.3$ keV is added to account for the iron line component. The fitting is further significantly improved with an additional hot thermal component ($kT\sim 4.1$ keV). The reduced chi-square decreases from $\chi^{2}/dof=1105.7/1042$ to $\chi^{2}/dof=1047.8/1040$, and the fit residuals become flat in the whole energy band (bottom panel of Figure 2). The same situation occurs for the NuSTAR data obtained at the outburst decay phase. Alternatively, the blackbody component is required with a confidence level of 98% according to $F$-test, but the iron line is too weak to be detected in the last NuSTAR observation. We, therefore, suggest that the NuSTAR spectra can be described as a combination of a hot blackbody and a cut-off power-law component, and the iron line emission is required at the high luminosity state (Table 2).

For the XMM–Newton data, we generate the spectral response files with the sas tasks rmfgen and arfgen, and rebin the spectra by using the task specgroup to have at least 20 counts per bin to enable the use of chi-square statistics and not to oversample the instrument energy resolution by more than a factor of three. La Palombara et al. (2018) carried out a detailed analysis on the XMM–Newton data, and concluded that the continuum spectrum was dominated by the power-law component, and displayed a soft excess below 2 keV. The latter feature was further described with the sum of a cool blackbody and a hot thermal plasma component. Here, we fit the pn and MOS1/2 data simultaneously and confirm that all three components are required by the data. Adopting the same model as used in La Palombara et al. (2018) [tbabs*(apec+bbodyrad+powerlaw+gaussian) in xspec] with the metal abundance of the APEC component fixed to 0.2 $Z_{\odot}$, we obtain the similar values for all parameters: $nH=0.07_{-0.02}^{+0.03}\times 10^{22}$ cm*-2*, $kT_{\rm apec}=0.96_{-0.09}^{+0.08}$ keV, Norm${}_{\rm apec}=6.4_{-2.1}^{+2.0}$, $kT_{\rm BB}^{\rm low}=0.22_{-0.02}^{+0.03}$ keV, Norm${}_{\rm BB}^{\rm low}=86.3_{-47.8}^{+90.4}$, $\Gamma=0.73_{-0.02}^{+0.02}$, Norm${}_{\rm PL}=2.66_{-0.09}^{+0.09}\times 10^{-3}$, $E_{\rm Gau}=6.34_{-0.25}^{+0.24}$ keV, $\sigma=0.34_{-0.21}^{+0.28}$ keV, Norm${}_{\rm Gau}=3.25_{-1.89}^{+2.31}\times 10^{-5}$, and $\chi^{2}/dof=575.2/432$). But because the peak emission of the hot blackbody component needed by the NuSTAR data is beyond 10 keV (Figure 3), its parameters cannot be constrained with the XMM–Newton spectra alone.

Since the separation of first NuSTAR observation and the XMM–Newton observation is less than 2 d, we also try to fit them together with the common model of tbabs*(apec+bbodyrad+bbodyrad+cutoffpl+gaussian) as discussed above. The multiplicative constant for pn data is frozen at unity, and those for MOS1/2 and FPMA/B are allowed to vary. The derived constant factors for MOS1/2, FPMA, and FPMB are $0.97\pm 0.01$, $1.15\pm 0.01$, and $1.22\pm 0.01$, respectively. The unfolded spectra are plotted in Figure 3, and the spectral parameters are presented in Table 2 and Figure 4. The hot thermal plasma, the two blackbody and the non-thermal components contribute the X-ray emissions (in 0.5–79 keV) of $\sim$ $2.4\times 10^{35}$ erg s*-1*, $1.3\times 10^{36}$ erg s*-1*, $6.9\times 10^{36}$ erg s*-1*, and $1.02\times 10^{38}$ erg s*-1*, respectively. There is no obvious evidence of cyclotron absorption line feature in either the XMM–Newton nor the NuSTAR data. Finally, we would caution that the small discrepancy between the joint-fitting results and those obtained from fitting the XMM–Newton and the NuSTAR data alone could be due to the calibration differences and the spectral evolution within 2 d.

In order to verify the existence of the blackbody components, we also try to use other commonly used models to describe the non-thermal X-ray continuum (e.g. Coburn et al., 2002; West et al., 2017), such as the negative and positive exponential cut-off (the so-called NPEX model, Mihara et al., 1998), the Fermi Dirac cut-off (the so-called FDCut model, Tanaka, 1986), and the high-energy cut-off power-law models (highecut*powerlaw in XSPEC). These models predict that the spectra having a power-law profile below 10 keV and rolling off in different ways at high-energy band. The soft excess revealed below 2 keV is not sensitive to the adopted non-thermal models. On the other hand, the temperature of hot thermal component does not change much while its emitting size could vary by a factor of $<3$ when different continuum models are used. That is much smaller than the radius of a NS. Alternatively, a more physical model, CompTT, is also used to fit the spectra resulting in the similar parameter values for the thermal component, but obtain a worse fit. Note that, compared to the cut-off power-law model, these models have more parameters, which sometimes are difficult to be constrained. In sum, we suggest that the spectral parameters yielded by the cut-off power-law model are reliable and can be better constrained.

2.3 Pulse profiles and pulse fractions

The 0.3–12 keV source events are extracted from XMM–Newton EPIC data and are barycentrically corrected with the command barycen. Meanwhile, for the NuSTAR data, the source events are extracted in 3–79 keV for the period calculation. For each observation, an accurate template profile with 50 phase bins is created by folding the whole event data. Then we divide one observation into several segments having equal exposure (4000 s), and derive the pulse times of arrivals (TOAs) of the pulsar by comparing the template profile with the one from each segment, as detailed in the following: (1) search for the best spin frequency using the Pearson $\chi^{2}$ method; (2) fold the pulse profile with the starting time of the observation as the reference epoch; (3) Calculate the phase shift using the cross-correlation between the pulse profile and the template profile, which represents the TOA of each observation. Finally, we determine the rotation frequencies and their derivatives for each observation by fitting the TOAs with TEMPO2 (Hobbs et al., 2006, Table 1).

We also produce the light curves in time resolution of 0.1 s. The barycentric corrected light curves are folded over the best-fitting period, and the pulse fractions are calculated as $PF=(C_{\rm max}-C_{\rm min})/(C_{\rm max}+C_{\rm min})$, where $C_{\rm max}$, $C_{\rm min}$ are the maximum and the minimum count rates of the profile. The evolutions of the pulse profile and pulse fraction are plotted in Figures 5 and 6, respectively. However, we cannot investigate the pulse modulation above 50 keV (30 keV) for the first two (the last) NuSTAR observations owing to the low count rate.

2.4 Results

Our results are summarized as follows: (I) During the 2017 giant outburst, SXP 59 reached a peak luminosity of $\sim 1.1\times 10^{38}$ erg s*-1*, that is 60% Eddington luminosity of a NS. (II) Investigating the XMM–Newton data, we confirm that the soft excess reported by La Palombara et al. (2018) consists of a cool thermal component ($kT_{\rm BB}\sim 0.2$ keV) with a size of $10^{2}$ km and a hot thermal plasma . (III) The hard X-ray spectra ($>3$ keV) are modelled by three components: a hot blackbody component ($kT_{\rm BB}\sim 1.5-4$ keV), a non-thermal component, and an iron emission line. The temperature of blackbody decreases with time, while its normalization remains constant ($\sim 0.01$, Figure 4), corresponding to a size of $R\sim 0.62$ km. (IV) The pulse profiles given by the first two NuSTAR data are energy dependent and have two narrow peaks at phase of $\sim$ 0.5 and 0.7. Alternatively, the pulse shape at the low luminosity state has a single peak profile. (V) The pulse fraction increases with the photon energy and saturates at 65% above 10 keV for all three NuSTAR observations.

3 Discussions AND Conclusions

The accretion geometry in BeXRBs is mainly governed by the NS magnetic field strength ($B$) and the accretion rate (Basko & Sunyaev, 1976; Riffert & Meszaros, 1988; Kraus et al., 1995; Becker et al., 2012; Mushtukov et al., 2018). For the case of low accretion rate, the falling material is funnelled by the magnetic field to small regions around the polar caps of NSs (i.e. hot spots). The X-ray flux is mainly contributed by the thermal component from the hot spots with a temperature of $>1$ keV and a small radius of $<1$ km, e.g. SAX J2103.5+4545 (İnam et al., 2004), 1A 0535+262 (Mukherjee & Paul, 2005), RX J1037.5-5647 (La Palombara et al., 2009). Theoretically, the size of the hot spot increases with luminosity (Lamb et al., 1973; Frank et al., 2002; Mushtukov et al., 2015). When the interaction between the thermal photons and the falling material (bulk motion Comptonization) is non-negligible, the observed spectrum would be deviated from the blackbody form, but can be fitted the CompTT model in xspec (e.g. Doroshenko et al., 2010; Tsygankov et al., 2019). As the accretion is larger than the critical value, the accretion column is formed and blocks the sight of hot spots. That is, the X-ray flux is dominated by the non-thermal component from the accretion column, and no emission from hot spots is expected.

The change of beam pattern (i.e. the existence of accretion column or hot spots at the stellar surface, the so-called fan beam and pencil beam) results in different pulse profiles. It has been observed in several giant outbursts of BeXRBs that, the pulse shapes transit from double peaks at high luminosity to single peak at low luminosity, and the pulse fraction increases with energy, e.g. 1A 0535+262 (Bildsten et al., 1997), SMC X-3 (Weng et al., 2017; Zhao et al., 2018), and Swift J0243.6+6124 (Tsygankov et al., 2018; Wilson-Hodge et al., 2018). Such evolution sequence can be interpreted as the different radiation beam patterns working in the super-critical and sub-critical accretion regimes (Basko & Sunyaev, 1976; Becker et al., 2012; Mushtukov et al., 2015). Taking account of the exact Compton scattering cross section in a strong magnetic field, Mushtukov et al. (2015) argued that the critical luminosity was not a monotonic function of $B$, and it reached a minimum of a few $10^{36}$ erg s*-1* when the cyclotron energy was about 10–20 keV (fig. 5 in their paper).

SXP 59 entered into a type II outburst in 2017 and became one of the brightest BeXRBs with a peak X-ray luminosity of $1.1\times 10^{38}$ erg s*-1*. Investigating the XMM–Newton and NuSTAR observations executed at different flux levels, we find that the pulse profiles evolve both with the photon energy and the X-ray luminosity (Figure 5). In general, the pulse profiles above 2 keV exhibit two narrow peaks at the high luminosity, and turn into a single peak in the last NuSTAR data. Although it is difficult to constrain changes in the geometry of emitting region with the data presented in this work, our results are in favor of the scenario that the source transited from the super-critical state to the sub-critical state as observed in 1A 0535+262 and SMC X-3. The critical luminosity is of $3.2\times 10^{36}~{}{\rm erg~{}s^{-1}}<L_{\rm crit}<5.9\times 10^{37}~{}{\rm erg~{}s^{-1}}$., that is a typical value for a Be X-ray pulsar (e.g. Becker et al., 2012; Mushtukov et al., 2015). It might further suggest a typical magnetic field ($\sim 10^{12}-10^{13}$ G) for the NS in SXP 59, although we cannot put tight constraint at current stage. It worth to note that the cyclotron absorption line feature is the direct evidence for the NS magnetic field; however, it could be transient and too weak to be detected. For instance, the bursting pulsar, GRO J1744-28 was discovered in 1995 (Fishman et al., 1995) and since then has been observed frequently by X-ray missions (e.g. BeppoSAX, RXTE, XMM–Newton, Chandra, and NuSTAR); but the weak absorption feature at $\sim$ 4.5 keV was detected only recently (D’Aì et al., 2015; Doroshenko et al., 2015). Therefore, the absence of cyclotron absorption line is not in contradiction with a typical magnetic field for SXP 59.

It was reported that, a cool thermal emission ($kT_{\rm BB}\sim 0.2$ keV) with a large emission area emerged in the XMM–Newton data of SXP 59 (La Palombara et al., 2018). The spectral modelling parameters along with the significantly small pulse fraction detected below 1 keV ($<35\%$, Figure 5), are in favor of the scenario that the central hard X-rays are reprocessed by the inner region of the accretion disc (Hickox et al., 2004; La Palombara et al., 2018). The X-ray continuum above 2 keV of SXP 59 is dominated by the non-thermal component from the accretion column, and can be phenomenologically fitted by a cut-off power-law component plus a hot blackbody emission. In this work, we do not find correlation between the parameters of the cut-off power-law model ($E_{\rm cut}$ and $\Gamma$) and the luminosity in the giant outburst of SXP 59 (Table 2 and Figure 4). On the other hand, the behavior of hot blackbody emission is quite puzzling. If the source is at the sub-critical state, this component is generally considered to be from the base of accretion column, and contributes a large portion of X-ray flux (e.g. La Palombara et al., 2009). However, the prediction that the hot spot shrinks by a factor of $\sim 3$ during the outburst decay of SXP 59, conflicts with the constant normalization derived from the data (Figure 4). On the other side, theoretically, we can not receive the hot spot emissions directly at the high luminosity due to the accretion column. Nevertheless, the hot blackbody component is needed to fit the spectra of some luminous Be X-ray pulsars ($L_{\rm X}>10^{36}-10^{39}$ erg s*-1*), e.g. GX 1+4 (Yoshida et al., 2017), EXO 2030+375 (Reig & Coe, 1999), SXP 59 (this work), in particular, Swift J2043.6+6124 (Tao et al., 2019). The unexpected hot blackbody emission challenges the canonical accretion theories, which are mostly based on a pure dipole magnetic field. These observational results would indicate that either more physical spectral models are required to describe the spectra of luminous X-ray pulsars, or that the magnetic filed configuration deviates from a dipole field close to the NS surface.

Acknowledgements

We thank the anonymous referee for the helpful comments, which significantly improve this work. SSW thanks Hua Feng and Lian Tao for valuable discussions. We acknowledge the use of public data from the High Energy Astrophysics Science Archive Research Center Online Service. This work is supported by the National Natural Science Foundation of China under grants 11673013, 11703014, 11503027, 11573023, U1838201, U1838104, and the Natural Science Foundation from Jiangsu Province of China (grant no. BK20171028).