Observing the transient pulsations of SMC X-1 with NuSTAR

TL;DR

This study uses NuSTAR observations to analyze transient pulsations in SMC X-1, revealing that pulsation variability is likely caused by changing obscuration rather than flux changes, and suggesting parallels with ultraluminous X-ray pulsars.

Contribution

It demonstrates that pulsation appearance in SMC X-1 is linked to obscuration effects, providing insights into the variability mechanisms of high-luminosity X-ray pulsars.

Findings

Pulsations switch on and off without large flux changes.

Spectral shape remains consistent across states, with variability due to obscuration.

SMC X-1's luminosity suggests it as an analog to ULXPs.

Abstract

We report on NuSTAR observations of transient pulsations in the neutron star X-ray binary SMC X-1. The transition from non-pulsing to pulsing states was not accompanied by a large change in flux. Instead, both pulsing and non-pulsing states were observed in a single observation during the low-flux super-orbital state. During the high-state, we measure a pulse period of $P = 0.70117(9)\,\mathrm{s}$ at $T_{ref} = 56145\,\mathrm{MJD}$. Spectral analysis during non-pulsing and pulsing states reveals that the observations can be consistently modeled by an absorbed power law with a phenomenological cutoff resembling a Fermi-Dirac distribution, or by a partially obscured cutoff power law. The shapes of the underlying continua show little variability between epochs, while the covering fraction and column density vary between super-orbital states. The strength of pulsations also varies, leading us to infer that the absence and reemergence of pulsations are related to changing obscuration, such as by a warped accretion disk. SMC X-1 is accreting near or above its Eddington limit, reaching an unabsorbed X-ray luminosity of $L_{\rm X}({\rm 2-10~keV}) \approx 5 \times 10^{38}\, {\rm erg}\, {\rm s}^{-1}$. This suggests that SMC X-1 may be a useful local analog to ultraluminous X-ray pulsars (ULXPs), which likewise exhibit strong variability in their pulsed fractions, as well as flux variability on similar timescales. In particular, the gradual pulse turn-on which has been observed in M82 X-2 is similar to the behavior we observe in SMC X-1. Thus we propose that pulse fraction variability of ULXPs may also be due to variable obscuration.

Figures

Tables

Table 1. Table 1: Results of the folding pulsation search for each epoch.

Epoch	$T_{ref}$ (MJD)	$P$ (s)	$\left|\dot{P}\right|\left({10}^{-8}s{s}^{-1}\right)$	Pulse Fraction (%)
I	$56113.28661210$	$\cdots$	$\cdots$
II	$56113.92279551$	$0.70121\left(20\right)$		$21.5\pm 1.5$
III	$56145.10372569$	$0.70117\left(9\right)$		$40.9\pm 0.5$

Table 2. Table 2: Values of spectral parameters determined by χ 2 superscript 𝜒 2 \chi^{2} fitting of observed spectra. Two models are shown: a fully covered power law with a Fermi-Dirac-like cutoff modeled by fdcut (top), and a partially covered power law with an exponential cutoff, modeled by cutoffpl (bottom).

Component	Parameter	Epoch I	Epoch II	Epoch III
tbabs	$N_H$ (${10}^{22}{\mathrm{cm}}^{-2}$)	$16\pm 5$	${24}_{-4}^{+5}$	${1.9}_{-0.9}^{+1.3}$
fdcut	$\Gamma$	$1.0$$\dagger$$\dagger$footnotemark:	$1.0$$\dagger$$\dagger$footnotemark:	$1.0$$\dagger$$\dagger$footnotemark:
	$E_{\mathrm{cut}}$ (keV)	${17.3}_{-2.3}^{+1.6}$	${11.0}_{-5.1}^{+3.2}$	${9.1}_{-3.0}^{+2.5}$
	$E_{\mathrm{fold}}$ (keV)	${6.7}_{-0.6}^{+0.8}$	${8.7}_{-0.8}^{+1.0}$	$9.6\pm 0.4$
	Norm (${10}^{-3}$)	${8.0}_{-0.9}^{+1.4}$	${13.8}_{-2.5}^{+3.2}$	${107}_{-16}^{+29}$
gauss	$E_{6.4}$ (keV)	$6.36\pm 0.04$	$6.36\pm 0.06$	${6.51}_{-0.07}^{+0.09}$
	${\sigma }_{6.4}$ (keV)	$0.24\pm 0.06$	$0.21\pm 0.10$	$0$$\dagger$$\dagger$footnotemark:
	Norm (${10}^{-4}$)	$3.3\pm 0.5$	${2.3}_{-0.5}^{+0.6}$	$1.9\pm 0.7$
gauss	$E_{13.5}$ (keV)	$13.5$$\dagger$$\dagger$footnotemark:	$13.5$$\dagger$$\dagger$footnotemark:	$\cdots$
	${\sigma }_{13.5}$ (keV)	${2.2}_{-0.7}^{+0.8}$	${1.7}_{-0.9}^{+1.8}$
	Norm (${10}^{-4}$)	${2.3}_{-1.1}^{+1.6}$	${1.1}_{-0.7}^{+1.7}$
bbody	$kT$ (keV)	${0.36}_{-0.06}^{+0.04}$	${0.31}_{-0.04}^{+0.03}$	$0.26\pm 0.09$
	Norm (${10}^{-3}$)	${3.6}_{-1.9}^{+3.5}$	${19.6}_{-9.3}^{+24.6}$	${56}_{-45}^{+195}$
bbody	$kT$ (keV)	$\cdots$	$\cdots$	$1.46\pm 0.07$
	Norm (${10}^{-3}$)			$2.3\pm 0.4$
	Absorbed Flux$a$$a$footnotemark: (${10}^{-11}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	${6.26}_{-0.32}^{+0.15}$	${6.67}_{-0.49}^{+0.06}$	${99.3}_{-15.4}^{+0.8}$
	Unabsorbed Flux$a$$a$footnotemark: (${10}^{-10}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	${1.46}_{-0.33}^{+0.56}$	${2.82}_{-0.69}^{+1.16}$	${11.6}_{-0.9}^{+4.0}$
	Unabsorbed Luminosity (${10}^{37}\mathrm{erg}{s}^{-1}$)	${6.43}_{-1.44}^{+2.48}$	${12.4}_{-3.0}^{+5.1}$	${51.0}_{-4.1}^{+17.7}$
	${\chi }^2/d.o.f.$	693/630 (1.10)	681/679 (1.00)	1091/1019 (1.07)
tbpcf	$N_H$ (${10}^{23}{\mathrm{cm}}^{-2}$)	${7.8}_{-2.1}^{+1.9}$	${5.6}_{-1.2}^{+0.9}$	${26}_{-10}^{+18}$
	$f_{\mathrm{covering}}$ (%)	${51}_{-8}^{+7}$	${61}_{-4}^{+2}$	${15}_{-6}^{+7}$
cutoffpl	$\Gamma$	$0.5$$\dagger$$\dagger$footnotemark:	$0.5$$\dagger$$\dagger$footnotemark:	$0.5$$\dagger$$\dagger$footnotemark:
	$E_{\mathrm{cut}}$ (keV)	$8.4\pm 0.6$	${8.9}_{-0.3}^{+0.4}$	$9.1\pm 0.2$
	Norm (${10}^{-3}$)	${6.7}_{-1.2}^{+1.5}$	${7.5}_{-0.9}^{+0.7}$	$53\pm 4$
gauss	$E_{6.4}$ (keV)	${6.35}_{-0.05}^{+0.04}$	${6.34}_{-0.07}^{+0.06}$	$6.52\pm 0.08$
	${\sigma }_{6.4}$ (keV)	$0.22\pm 0.07$	${0.24}_{-0.09}^{+0.10}$	$0$$\dagger$$\dagger$footnotemark:
	Norm (${10}^{-4}$)	${3.6}_{-0.6}^{+0.7}$	${2.6}_{-0.5}^{+0.7}$	${2.0}_{-0.9}^{+0.8}$
gauss	$E_{13.5}$ (keV)	$13.5$$\dagger$$\dagger$footnotemark:	$13.5$$\dagger$$\dagger$footnotemark:	$\cdots$
	${\sigma }_{13.5}$ (keV)	${4.9}_{-1.1}^{+1.3}$	${1.7}_{-1.1}^{+3.3}$
	Norm (${10}^{-4}$)	$7.8\pm 3.1$	${1.0}_{-0.7}^{+3.0}$
bbody	$kT$ (keV)	$\cdots$	$\cdots$	${0.23}_{-0.07}^{+0.08}$
	Norm (${10}^{-1}$)			${1.1}_{-1.0}^{+72.3}$
bbody	$kT$ (keV)	$\cdots$	$\cdots$	${1.46}_{-0.09}^{+0.08}$
	Norm (${10}^{-3}$)			${3.0}_{-0.7}^{+0.8}$
	Absorbed Flux$a$$a$footnotemark: (${10}^{-11}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	${6.49}_{-0.24}^{+0.06}$	${6.90}_{-0.09}^{+0.08}$	${105}_{-18}^{+2}$
	Unabsorbed Flux$a$$a$footnotemark: (${10}^{-10}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	${1.06}_{-0.14}^{+0.17}$	${1.18}_{-0.11}^{+0.08}$	${12.7}_{-0.6}^{+1.1}$
	Unabsorbed Luminosity (${10}^{37}\mathrm{erg}{s}^{-1}$)	${4.65}_{-0.62}^{+0.73}$	${5.20}_{-0.47}^{+0.34}$	${55.8}_{-2.8}^{+4.8}$
	${\chi }^2/d.o.f.$	705/632 (1.11)	681/681 (1.00)	1090/1019 (1.07)

Equations

$$PF=\frac{F_{\max }-F_{\min }}{F_{\max }+F_{\min }}$$

$$T_n=T_0+n{P}_{\mathrm{orb}}+\frac{1}{2}{n}^2{P}_{\mathrm{orb}}{\dot{P}}_{\mathrm{orb}}$$

$$f_{\mathrm{FD}}(E)=\frac{E^{-\Gamma }}{1+e^{(E-E_{\mathrm{cut}})/E_{\mathrm{fold}}}}$$

$$\mathrm{BIC}=k\ln (n)+{\chi }^2$$

Full text

accretion, accretion disks – stars: neutron – X-rays: binaries – X-rays: individual (SMC X-1)

††software: Stingray (Huppenkothen et al., 2016), HENDRICS (Bachetti, 2015), NUSTARDAS, MaLTPyNT (Bachetti, 2015), DS9 (Joye & Mandel, 2003)

1 Introduction

High Mass X-ray binary (HMXB) systems, in which a compact object accretes matter from a more massive companion, often exhibit variability on multiple timescales ranging from less than a second to several months. For most HMXBs, accretion onto the compact object is fed by the stellar winds of the companion. However, in some cases, the objects may orbit at a small enough distance that the companion fills its Roche lobe, resulting in higher accretion rates and higher luminosities. In the case of accreting neutron stars, matter is funnelled along the magnetic field lines onto the surface of the star, resulting in a column of accreted material at the magnetic poles. Because the magnetic poles and the spin axes are not perfectly aligned, the accretion column revolves at the same rate as the neutron star, resulting in emission which appears pulsed.

One such HMXB which has been well-studied is SMC X-1. This system, residing in the Small Magellanic Cloud, was first detected by Price et al. (1971). It was later resolved as a discrete source by Leong et al. (1971), who reported significant variability in both the intensity and spectrum of the source. The binary nature of SMC X-1 was soon confirmed by Schreier et al. (1972) who discovered periodic occultations with an orbital period of around 3.9 days. SMC X-1 also exhibits pulsations with a period of about 0.7 seconds, and the pulse fraction and shape are known to vary significantly over time (Lucke et al., 1976). The existence of X-ray pulsations confirms that the accreting compact object is a neutron star. Accretion onto the neutron star has been attributed to Roche lobe overflow (Hutchings et al., 1977; van Paradijs & Zuiderwijk, 1977) of its companion, Sk 160, which has been spectrographically classified as a B0 I supergiant (Webster et al., 1972). This classification places SMC X-1 in a subcategory of HMXBs known as supergiant X-ray binaries (SGXB). Finally, the source exhibits super-orbital variability on a timescale of 45 to 60 days, which has been attributed to obscuration by a precessing tilted accretion disk (Wojdowski et al., 1998). Throughout this paper, we assume a distance to SMC X-1 of 60.6 kpc as reported by Hilditch et al. (2005).

One of only a handful of SGXBs known to accrete via Roche lobe overflow, SMC X-1 exhibits persistent emission near or above its isotropic Eddington luminosity of $L_{\rm Edd}\sim 1.3\times 10^{38}\,\mathrm{erg\,s^{-1}}$ (for a mass of $\sim 1.1\,M_{\odot}$, as reported by van der Meer et al. 2007) varying between $L_{\rm X}(2-12\,{\rm keV})\sim 10^{37}\ \mathrm{erg}\ \mathrm{s}^{-1}$ in the low state and luminosities in excess of $5\times 10^{38}\ \mathrm{erg}\ \mathrm{s}^{-1}(2-12\,{\rm keV})$, more than three times its Eddington luminosity, in the high state (Bonnet-Bidaud & Klis, 1981). In addition to this persistent emission, SMC X-1 has been shown to exhibit type II X-ray bursts with durations of tens of seconds (Angelini et al., 1991; Rai et al., 2018). Its near- to super-Eddington luminosity places the source in a middle ground between less luminous Be/X-ray binaries (BeXB), which exhibit a range of persistent X-ray luminosities from $10^{32}\ \mathrm{erg}\ \mathrm{s}^{-1}$ (Tomsick et al., 2011) up to $10^{35}\ \mathrm{erg}\ \mathrm{s}^{-1}$ (Reig & Roche, 1999), and brighter ultraluminous X-ray pulsars (ULXPs).

ULXPs, the known examples of which are M82 X-2 (Bachetti et al., 2014), NGC 7793 P13 (Fürst et al., 2016; Israel et al., 2017a), NGC 300 ULX1 (also SN 2010da, Carpano et al. 2018) and NGC 5907 ULX-1 (Israel et al., 2017b), vary between bright pulsing states during which the luminosity can reach $10^{41}\ \mathrm{erg}\ \mathrm{s}^{-1}$ — several orders of magnitude above the Eddington limit of a typical neutron star — and faint states when the luminosity drops to $10^{37-38}\ \mathrm{erg}\ \mathrm{s}^{-1}$ (Kaaret et al., 2017). Similar to SMC X-1 (Inam et al., 2010), these sources exhibit pulsations with periods on the order of one second and spin-up rates of $|\dot{P}|=10^{-11}-10^{-9}\,\mathrm{s\,s^{-1}}$, with the exception of NGC 300 ULX1, which has a much longer pulse period of 32 s and a faster spin-up rate on the order of $10^{-7}\,\mathrm{s\,s^{-1}}$.

A given ULXP may not exhibit detectable pulsations at all times, and when they are detected, the fraction of their flux which is pulsed is variable. Pulse transience has been attributed to the propeller effect in which rotation of the neutron star’s magnetosphere halts accretion by flinging accreting material out of the system before it can reach the corotation radius (Illarionov & Sunyaev, 1975). In contrast to the flux variability of SMC X-1 which occurs quasi-periodically with a continuous transition between high and low states, the propeller effect is associated with changes of more than a factor of $40$ in luminosity on shorter timescales which results in a bimodal flux distribution in ULXPs (Tsygankov et al., 2016). In terms of pulse fraction, this bimodality corresponds to distinct pulsed and non-pulsed states. However, continuous variability in pulse fraction has also been observed in ULXPs. In particular, the pulse fraction of M82 X-2 was shown to gradually increase from $8$% to $23$% in the $10-30$ keV range over an interval of around 10 days (Bachetti et al., 2014).

Periodic variability on timescales of 60-80 days has also been measured for the ULXPs NGC 7793 P13, NGC 5907 ULX-1, and M82 X-2 (Motch et al., 2014; Walton et al., 2016; Brightman et al., 2019, respectively). While the  64 day period observed in NGC 7793 P13 has been attributed to the orbital motion of the binary (Fürst et al., 2018), this variability has been classified as super-orbital in the case of M82 X-2. It is still uncertain whether the 78 day period observed in NGC 5907 ULX-1 is orbital or super-orbital in nature. Given that SMC X-1 displays super-orbital modulations on similar timescales, as well as its persistent near- to super-Eddington luminosity and variable pulsations, SMC X-1 may provide a link between ULXPs and classes of X-ray binaries which have been studied in more detail.

In this paper, we present timing and spectral analyses of two observations of SMC X-1. In Section 2, we describe the observations of SMC X-1 and our data reduction methods, including data extraction and corrections. In Section 3, we describe the methods and results of our timing analysis, and in Section 4, we describe the spectral analysis of SMC X-1. In Section 5, we discuss the results of our analyses and offer a physical interpretation. Finally, in Section 6, we list our conclusions and discuss possible applications of our analysis to studies of ULXPs.

2 Observations and Data Reduction

SMC X-1  was observed twice by the Nuclear Spectroscopic Telescope Array (NuSTAR) (Harrison et al., 2013) in 2012, in the first two months after the launch of the satellite for the purpose of calibration. NuSTAR consists of two focal plane modules, FPMA and FPMB, each of which is made up of four pixelated detectors (DET0-DET3). Each module has a field of view of about 12 arcminutes, and, combined with focusing optics at a focal length of 10 m, achieves an angular resolution of 18 arcseconds, full width at half maximum (FWHM). The energy resolution, given by the FWHM, is 400 eV at 10 keV and 900 eV at 68 keV, and the full energy range is 3-79 keV. The timing resolution of the onboard clock is 2 $\mu$s with a dead time of 2.5 ms, leading to a maximum count rate of around 400 events $\mathrm{s^{-1}}$.

The first observation took place on 2012 July 5 (OBSID 10002013001) and the second took place on 2012 August 6 (OBSID 10002013003) with exposure times of 27 ks and 15 ks, respectively. Figure 1 shows the light curve of the source as observed by the Neil Gehrels Swift Observatory Burst Alert Telescope (BAT) during a 250 day interval bracketing the NuSTAR observations in 2012. The super-orbital period of around 60 days is clearly visible, and the red bars show the location and duration of each observation in the super-orbital cycle. The first NuSTAR observation occurred at the end of the low state, when the luminosity was just beginning to rise, while the second observation occurred near the end of the high state, when the source was growing fainter. The observations were planned such that they avoided obscuration effects due to the donor star.

We reduced the data using version 1.8.0 of the NuSTARDAS pipeline and NuSTAR CALDB v20170817. We used DS9 (Joye & Mandel, 2003) to select a circular source region with radius 55 arcseconds centered on the position of the source determined by automatic centroid detection. We also selected a circular background region with radius 80 arcseconds located on the same detector as the source, taking care to choose a region free of other sources and outside the source distribution. We corrected the photon arrival times to the solar system barycenter using the position of the source used for data extraction. Before analysis, the photon arrival times were also corrected for the orbital motion of the source using parameters reported by Falanga et al. (2015) and Inam et al. (2010).

We define three epochs of observation, labeled Epochs I, II, and III. The NuSTAR light curve for each observation is shown in Figure 2. Epoch I is defined as the first 40 ks (13 ks of exposure time) of observation 10002013001, while the latter half (14 ks of exposure) of observation 10002013001 makes up Epoch II. The whole of observation 10002013003 makes up Epoch III, which has an exposure time of 15 ks. During observation 10002013001, the source was positioned on DET0, while the source was positioned on DET3 near the gap between DET3 and DET0 during observation 10002013003. Movement of the source between the two detectors accounts for the $\sim 5000$ s variability apparent in Figure 2b. The background count rate did not vary significantly between observations, and for all three epochs, the background rate remained below $10\%$ of the total count rate for energies up to $\sim 50$ keV. To avoid background contamination, we performed spectral analysis for energies between 3 keV and 40 keV, resulting in $5.2\times 10^{4}$, $6.3\times 10^{4}$, and $6.3\times 10^{5}$ spectral counts (combined FPMA and FPMB) for Epochs I, II, and III, respectively. For the purpose of spectral analysis, we binned the data such that there are at least 50 counts in each energy bin in Epoch I and Epoch II, and at least 100 counts in each energy bin in Epoch III. We chose to bin Epoch III with more events per bin due to the significantly higher count rate during that epoch.

3 Timing Analysis

We performed a timing analysis of both observations using the Stingray (Huppenkothen et al., 2016) and HENDRICS (Bachetti, 2015) software packages in order to determine the pulse fraction, pulse period, and spin-up rate during each epoch. The results of this analysis are shown in Table 1. The pulse fraction, $PF$, is defined as follows

[TABLE]

where $F_{\rm max}$ and $F_{\rm min}$ are the maximum and minimum fluxes in the pulse profile, respectively. All pulse fractions and corresponding errors quoted were calculated using a Monte Carlo analysis. Given a measured pulse period and derivative, we folded the observed events into a pulse profile with sixteen phase bins per cycle. The uncertainty in flux for each phase bin is given by a Poisson distribution. We sampled this distribution for each phase bin to produce a large number of simulated pulse profiles and passed these profiles through a Savitzky-Golay filter (Savitzky & Golay, 1964). We thus arrived at a distribution of smoothed profiles from which we extracted the mean pulse fraction and corresponding confidence regions. All uncertainties and upper limits quoted in this section and following sections correspond to 90% confidence ranges unless otherwise indicated.

Before searching for and analyzing pulsations, we first determined the orbital phase of the observations. Using the orbital parameters reported by Falanga et al. (2015), we determined the mid-eclipse times which occurred immediately before and after the observations according to the quadratic orbital change function

[TABLE]

where $T_{0}$ is the reference epoch (MJD 52846.6888), n is the number of elapsed orbits, $P_{\rm orb}$ is the orbital period measured at $T_{0}$ (3.8919232 days), and $\dot{P}_{\rm orb}$ is the time derivative of the orbital period measured at $T_{0}$ ($\mathrm{-3.77\times 10^{-8}\,day\,day^{-1}}$). Defining the mid-eclipse times preceding and following each observation as orbital phases 0 and 1, respectively, we found that the first observation occurred between orbital phases $0.34$ and $0.61$, while the second observation occurred between phases $0.52$ and $0.61$. These orbital phases are determined to better than $10^{-5}$ and lie far from the eclipse ingress and egress times. Therefore we can be confident that there were no obscuration effects due to the supergiant companion.

We performed pulsation searches on the combined filtered and calibrated FPMA and FPMB events for each of the three epochs. When combining the FPMA and FPMB events for each epoch, we produced common good-time intervals (GTIs) in order to avoid introducing artificial variability due to non-simultaneous observation and differences in sensitivity between the two focal plane modules. We began our pulsation search by performing a dynamic search using the HENDRICS function dyn_folding_search. This function steps over time and pulse period, folds the events into a pulse profile with that period, and calculates the $Z_{4}^{2}$ statistic (Buccheri et al., 1983), a measure of the probability of pulsation detection, of the profile produced at each step. The probability density function of the $Z_{4}^{2}$ statistic is equivalent to that of a $\chi^{2}$ distribution with 8 degrees of freedom. Therefore, one can use the $\chi^{2}$ cumulative distribution function to determine the probability that a pulse profile with a given value of $Z_{4}^{2}$ has been produced by noise. For example, a profile with $Z_{4}^{2}=13$ has a 10% probability of being produced by noise, therefore this can be considered a detection at 90% confidence. A 5-sigma detection, corresponding to a probability of $5.7\times 10^{-7}$ that a signal has been produced by noise, would yield a $Z_{4}^{2}$ statistic of 44.

The results of the dynamic pulsation search applied to Epoch III (Figure 3g) confirm the presence of strong pulsations. The pulsations appear to remain persistent throughout the observation with a period around $0.701\,\mathrm{s}$ and without a large period derivative. The results of this test are less striking upon application to Epochs I and II (Figures 3a and 3d). The $Z_{4}^{2}$ statistic reaches only a fraction of the maximum value measured during Epoch III, and during Epoch I, there is no sign of pulsations. However, during Epoch II, pulsations appear to have begun with a period similar to that observed in Epoch III, reaching a maximum detection probability at the end of the observation.

We next simultaneously searched for the period and first period derivative of the pulsations for each of the three epochs using the HENDRICS function folding_search. The results of this search are shown in the second column of Figure 3. We were able to measure the pulse period and to put upper limits on the first derivative during Epochs II and III by fitting the resulting $Z_{4}^{2}$ distributions to 2-dimensional Gaussian distributions111The uncertainties reported for the pulse periods and the upper limits of the spin-up rates were determined using the widths of the fitted Gaussian distributions. with mean pulse periods at $P_{\mathrm{II}}=0.70121(20)\,\mathrm{s}$ and $P_{\mathrm{III}}=0.70117(9)\,\mathrm{s}$, respectively. We have obtained upper limits on the instantaneous spin-up rates of $|\dot{P}_{\mathrm{II}}|<1.2\times 10^{-8}\,\mathrm{s\,s^{-1}}$ and $|\dot{P}_{\mathrm{III}}|<7.7\times 10^{-9}\,\mathrm{s\,s^{-1}}$. Note that there is a correlation between $P$ and $\dot{P}$ apparent in Figure 3. This is not intrinsic to the source itself but is an artifact introduced by the search procedure. In addition to measuring the pulse periods and constraining the instantaneous spin-up rates, we have also placed an upper limit on the secular spin-up rate between Epoch II and Epoch III of $|\dot{P}_{\rm sec}|<10^{-10}\mathrm{\,s\,s^{-1}}$. After determining the pulse periods and the spin-up rates, we then folded the events into pulse profiles at the $Z_{4}^{2}$ maxima produced by the pulsation searches. These pulse profiles are shown in Figures 3f and 3i. We observe distinct pulsations in the pulse profiles of Epochs II and III. The probability that these profiles were produced by noise is vanishingly small, being less than $10^{-37}$ in both cases.

In stark contrast, we were completely unable to constrain the pulse period during Epoch I. There are multiple local maxima of comparable amplitude in the $Z_{4}^{2}$ distribution. We therefore chose to fold the events into a pulse profile (see Figure 3c) using the maximum nearest the values measured during Epochs II and III. This corresponds to a pulse period of $P_{\rm fold}=0.70113161\,\mathrm{s}$ and a first derivative of $\dot{P}_{\rm fold}=2.69\times 10^{-9}\,\mathrm{s\,s^{-1}}$. The resulting profile has a pulse fraction of $<4.5\%$. This is relatively small compared to the pulse fractions of $21.5\%\pm 1.5\%$ during Epoch II and $40.9\%\pm 0.5\%$ during Epoch III. In addition to the small pulse fraction during Epoch I, the $Z_{4}^{2}$ value of the calculated pulse profile is less than 15 and corresponds to a probability of 7% that the detection is due to noise. Furthermore, when the last 5000 s of Epoch I are omitted from the pulsation search even this weak detection disappears, indicating that pulsations were absent until the very end of Epoch I. Therefore, we refer to Epoch I as the non-pulsing state.

The pulse periods that we have measured during Epochs II and III and the resulting pulse profiles are in line with previous measurements (cf. Moon et al., 2003; Naik & Paul, 2004; Raichur & Paul, 2010; Inam et al., 2010). In particular, we have extrapolated previous results by applying an orthogonal distance regression to the pulse frequencies reported by Inam et al.. We arrived at a spin up of $\dot{f}_{pulse}=2.589(8)\times 10^{-11}\mathrm{\,Hz\,s}^{-1}$ during the interval 50093-52988 MJD. When propagated forward to the beginning of Epoch III, a pulse period of $0.70093(2)\,\mathrm{s}$ is predicted. The discrepancy of $2.39(96)\times 10^{-4}$ s is small but nonzero. This is consistent with a piece-wise spin-up evolution, reported by Inam et al., in which the spin-up rate is variable. We also note that, although the pulse fraction increases with energy, the shapes of the pulse profiles during Epochs II and III do not appear to vary significantly with energy.

4 Spectral Analysis

We also performed a spectral analysis of each of the three epochs using Xspec (v.12.10.0 Arnaud, 1996). We simultaneously fit the spectra measured by FPMA and FPMB while including a relative constant to account for small ($<10\%$) differences in flux between the two focal plane modules. In addition, for all spectral models described in this section, we have included an absorber in the form of the tbabs component. This component compensates for absorption due to Galactic material. We fixed the equivalent HI column density of this component at $N_{\rm H}=4.58\times 10^{21}\ \mathrm{cm}^{-2}$, determined from the full-sky HI survey, HI4PI (Ben Bekhti et al., 2016). The spectral fits were performed using interstellar medium abundances reported by Wilms et al. (2000).

The spectra for each epoch are shown in Figure 4, and the results of our spectral analysis are presented in Table 2. Each panel in Figure 4 also includes residuals for three different models, including a simple absorbed power law meant to illustrate additional structure in the spectra. We have found two models that provide fits of similar quality and which result in physically reasonable parameters. Motivated by previous work by, e.g., Woo et al. (1995a), Angelini et al. (1991), and Pottschmidt et al. (2014), the first model we investigated was an absorbed power law with a phenomenological cutoff, named fdcut (Tanaka, 1986) for its resemblance to the Fermi-Dirac distribution, which has both a cutoff energy and folding energy and can be written

[TABLE]

where $\Gamma$ is the photon index, $E_{\rm cut}$ is the cutoff energy, and $E_{\rm fold}$ is the folding energy. The absorber in this model is fully covering and is modeled by tbabs (Wilms et al., 2000). The second model consists of a power law with an exponential cutoff, represented by the Xspec model cutoffpl, partially covered by an absorber modeled by tbpcf. In addition to these base models, we found that the fits benefited from the addition of secondary components, differing depending on the epoch. Below, we describe each of these models in more detail.

In order to compare the usefulness of additional model components, we use the Bayesian Information Criterion (BIC; Schwarz, 1978). In the case of $\chi^{2}$ fitting in Xspec, the BIC is given by

[TABLE]

where $n$ is the number of PHA bins being fitted and $k$ is the number of parameters estimated by a given model. For a given data set, model selection can be achieved by minimizing the BIC, which penalizes models with many parameters. For our analysis, $n$ lies between 600 and 1100 bins, meaning that removing one parameter from a model without a change in $\chi^{2}$ results in a decrease in the BIC of $\Delta\mathrm{BIC}\approx-7$. In determining the impact of adding or subtracting components, this may be considered one “unit” of model improvement.

The spectra observed during the two observations are qualitatively different, as is visible in Figure 4. For Epochs I and II, an absorbed fdcut model alone results in significant excess residuals below 4 keV, around 6.4 keV, and above 10 keV. The excess around 6.4 keV is consistent with previous detections of an Fe K$\alpha$ line in SMC X-1, such as those by Woo et al. (1995b) and Naik & Paul (2004). We included a Gaussian component at this energy to model the line, allowing both the position and width of the line to vary. To address the low-energy excess, we added a black body component with temperature $kT_{\rm BB}<0.5$ keV. Such a component has previously been detected in observations of SMC X-1 by the Chandra X-ray Observatory and XMM-Newton (Neilsen et al., 2004; Hickox & Vrtilek, 2005). Each of these components decreases the BIC by about 100, and adding both of these components results in a combined improvement to the fit of $\Delta\mathrm{BIC}=-304$ for Epoch I and $\Delta\mathrm{BIC}=-214$ for Epoch II, indicating that the improvement is significant. Adding these components does not resolve the “bump” above 10 keV. This residual resembles the 10 keV feature observed in other accreting pulsars (e.g. Coburn et al., 2002; Mihara, 1995; Santangelo et al., 1998), leading us to include a Gaussian component at $E=13.5$ keV. We froze the position of this component in order to better constrain other parameters, while the width of the Gaussian was allowed to vary. Adding this component does not result in a significant improvement to the fit, with $\Delta\mathrm{BIC}=-7$ for Epoch I and $\Delta\mathrm{BIC}=+6$ for Epoch II. However, because the residuals are clearly reduced, and this feature is consistent with previous studies of accreting pulsars, the component is included in the final fit.

In the case of Epoch III, the fdcut model again requires a low-temperature ($kT_{\rm BB}<0.5$ keV) black body in order to explain excess flux below 4 keV, the addition of which improves the fit by $\Delta\mathrm{BIC}=-93$. In addition, although the excess near 6.4 keV is not as prominent as in the previous two Epochs, adding a line at this energy improves the fit significantly ($\Delta\mathrm{BIC}=-46$). However, the width of the line is poorly constrained, leading us to freeze it at $\sigma=0$. We also found that for Epoch III, adding a $kT\approx 1.5$ keV blackbody component, like the one included by Pottschmidt et al. (2014) in their analysis of this observation, improves the fit by $\Delta\mathrm{BIC}=-46$ while also eliminating an excess of flux above 20 keV. Similar blackbody components with temperatures ranging between 1.2 keV and 3 keV have proved useful for modeling the spectra of several BeXRBs (e.g. Reig & Roche, 1999; La Palombara et al., 2009; Caballero et al., 2013). In contrast to Epochs I and II, the 13.5 keV bump is not observed during Epoch III, and its addition to the fdcut model does not improve the fit nor can this component be easily constrained.

The partially covered cutoff power law provides a similarly good fit to the data as the Fermi-Dirac-like model. However, the secondary components differ somewhat. Adding a line at $E=6.4$ keV to the base model again improves the fit significantly during Epochs I and II, with $\Delta\mathrm{BIC}=-97$ and $\Delta\mathrm{BIC}=-43$, respectively. Adding this line to Epoch III does not result in a striking improvement (only $\Delta\mathrm{BIC}=-6$), but the position of the Gaussian is constrained to the same value as in the fdcut model. The 13.5 keV bump is also added to this model for Epochs I and II, again slightly improving the fits. Unlike the absorbed fdcut model, the partially covered cutoff power law does not require a low-temperature black body to resolve excess emission below 4 keV during Epochs I and II. This component, along with the $kT_{\rm BB}\approx 1.5$ keV black body, remains in Epoch III. Adding each of these black body components individually yields different results. Including only the warm $kT_{\rm BB}\approx 1.5$ keV black body improves the fit by $\Delta\mathrm{BIC}=-142$, while adding only the low-temperature $kT_{\rm BB}<0.5$ keV black body does not improve the fit, yielding $\Delta\mathrm{BIC}=+6$. However, when both components are included, the fit is improved by $\Delta\mathrm{BIC}=-170$. In other words the combination of the low-temperature and high-temperature black body components improves the fit more than each of these components individually.

None of the three spectra can be fit to a simple one-component model, instead requiring several secondary components in order to properly fit the NuSTAR observations. In order to reduce degeneracies resulting from the number of parameters used in the final models, we froze some key parameters at values which are consistent with initial estimates. As mentioned above, the position of the bump above 10 keV was frozen at $E=13.5$ keV, and the width of the Fe K$\alpha$-like line in Epoch III was frozen at $\sigma=0$ keV. In addition to these, we froze the photon indices across all three Epochs at $\Gamma=1.0$ for the fdcut model and $\Gamma=0.5$ for the cutoffpl model.

The final model components, parameter estimates, and fit information are shown in Table 2. Here we remind the reader that the uncertainties quoted on spectral parameters represent 90% confidence regions. We found that for both models, the absorption parameters vary between epochs, while the underlying power law models show less variability. In the case of the fdcut model, the absorbing column density decreases by an order of magnitude from $N_{\rm H}=(2.4^{+0.5}_{-0.4})\times 10^{23}\,\mathrm{cm}^{-2}$ during Epoch II, to $N_{\rm H}=(1.9^{+1.3}_{-0.9})\times 10^{22}\,\mathrm{cm}^{-2}$ during Epoch III. As we have shown, the pulse fraction simultaneously increases between these two epochs. The shape of the fdcut component on the other hand remains consistent between Epochs II and III, with the cutoff energy of $E_{\rm cut}\approx 10\,{\rm keV}$ folding energy of $E_{\rm fold}\approx 9\,{\rm keV}$. However, Epoch I has a slightly higher cutoff energy and lower folding energy: $E_{\rm cut}=17.3^{+1.6}_{-2.3}\,{\rm keV}$ and $E_{\rm fold}=6.7^{+0.8}_{-0.6}\,{\rm keV}$.

On the other hand, in the case of the partially covered cutoff power law, the shape of the cutoffpl component stays constant. The exponential cutoff is consistent with $E_{\rm cut}\approx 9\,{\rm keV}$ during all three epochs. The absorption parameters show little variation between Epoch I and Epoch II, but the covering fraction drops by a factor of four from $f_{\rm covering}=(61^{+2}_{-4})\%$ in Epoch II to $f_{\rm covering}=(15^{+7}_{-6})\%$ in Epoch III. Between these epochs, the column density appears to increase by a factor of a few, but this parameter is poorly constrained during Epoch III due to the low covering fraction.

In both models, the underlying continuum increases in flux between successive epochs while the flux of the apparent Fe K$\alpha$ line remains constant. Taken together, these observations indicate that the increase in total flux between epochs cannot be attributed solely to the absorption included in the models described above, and that the source of the Fe K$\alpha$ line is likely distinct from the source of the continuum (e.g., originating in the photoionization region surrounding the central X-ray source). In addition, the appearance of the $kT_{\rm BB}\approx 1.5\,{\rm keV}$ black body in Epoch III, observed in both the fdcut model and the cutoffpl model, may indicate that the emitting region responsible for this component either did not exist or was obscured during Epochs I and II.

5 Discussion

Our timing analysis has shown that the source was observed in a non-pulsing state during Epoch I which subsequently evolved into a pulsing state, observed in Epoch II. During Epoch III, about a month after Epoch II, the pulsations had increased in strength, with the pulse fraction increasing by nearly a factor of two. At the same time, our spectral analysis has shown that for all three epochs, the emission of the source can be described by two different models: a fully covered power law with a phenomenological Fermi-Dirac-like cutoff, and a partially absorbed power law with an exponential cutoff. Each of these models requires additional components, but we found that both models are consistent with variable absorption parameters between the low and high pulse fraction states. In particular, the Fermi-Dirac model exhibits a decrease in absorption column density between Epoch II and Epoch III, and the cutoff power law is consistent with a decrease in the covering fraction (and a poorly constrained increase in column density), between the low and high pulse fraction states. In addition, the luminosity was observed to gradually increase between the non-pulsing Epoch I and the pulsing Epoch II.

In order to synthesize these results, we propose that the pulsing region was observed emerging from behind absorbing material. Given that Epochs I and II took place near the end of the low state, as illustrated in Figure 1 and that the super-orbital period has been attributed to a warped precessing accretion disk (Wojdowski et al., 1998; Clarkson et al., 2003; Dage et al., 2018), the absorbing material obscuring the pulsing region is likely part of the accretion disk. In short, the warped accretion disk absorbs and scatters the pulsed emission from the neutron star, leading to the absence of detected pulsations in the low state and a gradual turn-on of pulsations as the disk moves out of the line of sight. This picture is consistent with the opaque inner disk region described by Hickox & Vrtilek (2005) to explain apparent reprocessing of pulsed emission. Their analysis describes the case when both the neutron star and the inner regions of the warped accretion disk are visible to the observer, while ours describes the opposite case when the warped disk lies between the neutron star and the observer, obscuring the pulsing emission regions.

The relatively high absorption column density of the partially covered cutoff power law in Epoch III does not immediately fit within this interpretation. Although this column density is not particularly well constrained, it is still well above the values measured during the first two epochs. This increased column density is accompanied by an increase in the brightness of the power law itself; in other words, the increased flux between the two observations cannot be attributed solely to the absorption included in the models presented here. Thus one interpretation of the combination of a relatively high column density and a relatively low covering fraction is that much of the absorber is completely Compton thick during Epochs I and II. The partially covering absorber represented by tbpcf, then, only models the optically thinner regions of the accretion disk leading the covering fraction and column density to be underestimated for the first two epochs. During Epoch III, according to this interpretation, the source is observed through an overall less opaque region of the accretion disk so that a higher column density is measurable.

Pulse fraction variability, including pulse drop-out, has been observed in several other accreting pulsars. In some cases, this variability has been attributed to changes in accretion via the propeller effect (Illarionov & Sunyaev, 1975). These include HMXBs Vela X-1 and GX 301$-$2, which have been shown to exhibit off-states during which the sources drop in luminosity and pulsations are no longer detected (Kreykenbohm et al., 2008; Fürst et al., 2011). LMC X-4, in which pulse drop-out and turn-on have been observed during the high state (Brumback et al., 2018), presents a different case. Still others, such as the low mass X-ray binary Her X-1, exhibit pulse fraction variability attributed to obscuration by warped accretion disks (Kuster et al., 2001). Of these examples, the case of variable obscuration in Her X-1 is most analogous to the behavior we have observed in SMC X-1.

6 Conclusions

We have performed spectral and timing analyses of the accreting neutron star binary SMC X-1 for three separate epochs occurring during two NuSTAR observations. Our timing analysis confirmed that the source was observed in the midst of a turn-on of pulsations, which subsequently increased in strength before strong pulsations were observed a month later. Our spectral analysis, which showed variable absorber parameters and luminosity, led us to conclude that the non-pulsing state was due to obscuration of the pulsing region by a warped accretion disk, and that the gradual turn-on was due to the emergence of the pulsing emission from behind the disk.

Similarly to SMC X-1, ULXPs are also known to exhibit variability in their luminosities and pulse fractions. In particular, the gradual change in pulse fraction observed in the beginning of the 2014 observation of M82 X-2 (Bachetti et al., 2014) may share the same physical origin as the pulse fraction variability we have observed in SMC X-1. In that case, the super-orbital periods observed in ULXPs may be attributable to precessing accretion disks which periodically obscure the pulsing source, resulting in variability in the observed pulse fractions. Spectral and timing analyses at different points in the super-orbital cycles of known ULXPs, like the analysis we have carried out for SMC X-1, may help to illuminate the accretion mechanism and causes of variability in this recently discovered class of X-ray binary.

This work was supported under NASA grant No. NNG08FD60C, and made use of data obtained with NuSTAR, a project led by Caltech, funded by NASA and managed by NASA/JPL and has utilized the NUSTARDAS software package, jointly developed by the ASDC (Italy) and Caltech (USA). We would also like to thank the anonymous referee for providing helpful comments that helped to improve the quality of this paper.