NuSTAR DISCOVERY OF AN UNUSUALLY STEADY LONG-TERM SPIN-UP OF THE Be BINARY 2RXP J130159.6–635806

2RXP J130159.6–635806 regularly fell into the FOV of various X-ray telescopes thanks to extensive observational campaigns dedicated to PSR B1259–63 which is located only $9\buildrel{\,\prime}\over{.} 55$ away. This allows us to investigate the long-term evolution of the pulse period. We analyzed the XMM-Newton (Jansen et al. 2001) archival data from 2007–2011 using the procedures described by Chernyakova et al. (2005) and Science Analysis Software version 14.0.0. The list of selected XMM-Newton observations with corresponding period measurements are shown in Table 2.

Figure 2 presents the long-term evolution of the period, showing that the 2RXP J130159.6–635806 neutron star has undergone very strong and steady spin-up during the last 20 years. The first available value of spin period, measured in 1994, is 735 s (Chernyakova et al. 2005). The most recent measurements by NuSTAR, from 2014, show the period to be around 643 s (Table 1). This means that during the last ∼20 years the spin period decreased by ∼92 s, corresponding to a mean spin-up rate of $\sim 1.4\times {10}^{-7}$ s s⁻¹. Figure 2 shows a change in the average spin-up rate, first reported by Chernyakova et al. (2005).

We approximated the long-term period evolution with a linear function with one change in slope. A fit shows that the break occurred at MJD 51300 ± 217 (mid 1999) with a spin-up rate before and after the break of $(4.3\pm 2.7)\times {10}^{-8}$ s s⁻¹ and $(1.774\pm 0.003)\times {10}^{-7}$ s s⁻¹, respectively. As seen from the fit parameters, the spin-up rate becomes significantly higher after the break. As noted above, the inset of Figure 2 shows that the NuSTAR data points fit the long-term spin-up rate with remarkably high precision.

Similar behavior was observed for the X-ray pulsar GX 1 + 4, which showed steady spin up for more than a decade (see, e.g., Bildsten et al. 1997; González-Galán et al. 2012). However, GX 1 + 4 belongs to the subclass of accreting X-ray pulsars known as symbiotic X-ray binaries. For BeXRPs, persistent sources typically demonstrate pulse periods that are relatively stable. Examples include the population of Be systems in the Small Magellanic Cloud (Klus et al. 2014) and the well-known low luminosity Galactic system X Persei (Lutovinov et al. 2012). Transient BeXRPs show strong spin-up during Type I and Type II outbursts (Bildsten et al. 1997) with significant spin-down episodes in between (see, e.g., Postnov et al. 2015). Therefore, 2RXP J130159.6–635806 is a unique source among the BeXRPs because it demonstrates steady and high long-term spin-up with a relatively low and stable luminosity.

Pulsar pulse profiles and their evolution with luminosity and energy band depend on the geometrical and physical properties of the emitting regions in the vicinity of the neutron star. In Figure 3, the NuSTAR pulse profiles of 2RXP J130159.6–635806 are shown in three different energy bands: 3–10, 10–20, and 20–40 keV. The lower panel shows "soft" ((10–20)/(3–10) keV) and "hard" ((20–40)/(10–20) keV) hardness ratios.

Zoom In Zoom Out Reset image size

At all energies, the pulse profile can roughly be divided into one main peak at phases 0.0–0.5 and two smaller peaks at phases 0.5–0.75 and 0.75–1.0. The main feature that changes with energy is the depth of the minimum at phase 0.0. As shown in the lower panel of Figure 3, while the "hard" hardness ratio is almost constant, the "soft" hardness ratio shows a maximum at phase 0.0 due to an increase in the depth of the minimum at 3–10 keV. Such behavior is caused by differences in the source spectrum with pulse phase (see Section 4).

Figure 4 shows the pulsed fraction as a function of energy. The pulsed fraction is defined as $\mathrm{PF}=({I}_{\mathrm{max}}-{I}_{\mathrm{min}})/({I}_{\mathrm{max}}+{I}_{\mathrm{min}})$, where I_max and I_min are the maximum and minimum intensities in the pulse profile, respectively. Defined in this way, the pulsed fraction is very high (about 80%). The alternative way to characterize the pulsed fraction is the relative root mean square (rms), which can be calculated using the following equation:

$$\begin{eqnarray}&&\mathrm{rms}=\displaystyle \frac{{(\displaystyle \frac{1}{N}{\displaystyle \sum }_{i=1}^{N}{({P}_{i}-\langle P\rangle )}^{2})}^{\frac{1}{2}}}{\langle P\rangle },\end{eqnarray} \tag{ 1 }$$

where P_i is the background-corrected count rate in a given bin of the pulse profile, $\langle P\rangle$ is the count rate averaged over the pulse period, and N is the total number of phase bins in the profile (N = 30 in our analysis). The rms deviation obtained in this way reflects the variability of the source pulse profile in a manner that is not sensitive to outliers like the narrow features seen in the profile of 2RXP J130159.6–635806. Therefore, this quantity has a value of around 30%, which is much lower than the classically determined pulsed fraction and also independent of the energy band (see Figure 4).

Zoom In Zoom Out Reset image size

It is interesting to note that in contrast to the majority of X-ray pulsars (Lutovinov & Tsygankov 2009), 2RXP J130159.6–635806 does not show an increase in the pulsed fraction at higher energies. Such uncharacteristic behavior was previously observed for another persistent BeXRP—RX J0440.9+4431 (Tsygankov et al. 2012). On the other hand, Figure 3 shows that the pulsed fraction increases somewhat with energy if I_min is defined from the the pulse plateau rather than from the pulse minimum. In other words, the peak-to-plateau difference slightly grows with energy.

The observed 20-year strong and steady spin-up of 2RXP J130159.6–635806 reveals the existence of a long-term accelerating torque, which indicates that the binary interaction must lead to regular accretion over a decade-long timescale, although this could certainly be episodic (e.g., at periastron passages). The torque can be transferred by matter accreted from either the disk around the neutron star or a stellar wind from the optical counterpart. Unfortunately, there is no strong observational evidence allowing us to distinguish between these two different accretion channels. In both scenarios, this process is defined mainly by the mass accretion rate and the magnetic field strength (see, e.g., Ghosh & Lamb 1979).

Due to the unknown distance of 2RXP J130159.6–635806, the luminosity and mass accretion rate are highly uncertain; the magnetic field is also unknown since no cyclotron line is found in the energy spectrum (Section 4). However, some qualitative conclusions about the interaction between the accretion disk and the neutron star magnetosphere can be made from the noise power spectrum of the X-ray pulsar.

According to the "perturbation propagation" model, stochastic variations of viscous stresses in the accretion disk cause variations of the mass accretion rate (Lyubarskii 1997; Churazov et al. 2001). This, in turn, results in a specific shape of the Power Density Spectrum (PDS) of the emerging light curve. Namely, it will appear as a power law with slope –1 to $-1.5$ (but the exact value is not well established for X-ray pulsars) up to the break frequency, which is the highest frequency that can be generated in the accretion disk (Lyubarskii 1997).

In the case of a highly magnetized neutron star, this maximal frequency is limited by the Keplerian frequency at the magnetospheric radius, above which one can expect a cutoff in the source PDS (Revnivtsev et al. 2009). If the source stays in spin equilibrium (corotation), the cutoff frequency will coincide with the spin frequency of the pulsar, whereas, in the case of spin-up (increased mass accretion rate), the magnetosphere will be squeezed and additional noise will be generated at higher frequencies. If the mass accretion rate is known (i.e., the distance to the source is known), this property of the PDS can be used to estimate the magnetic field strength of the neutron star (Revnivtsev et al. 2009; Tsygankov et al. 2012; Doroshenko et al. 2014). The appearance of the break in the PDS does not necessarily indicate that accretion is from a disk. There are wind-accreting sources in spin equilibrium also showing a break in their PDSs around the pulse frequency (Hoshino & Takeshima 1993). However, the evolution of the PDS shape as a function of mass accretion rate in such systems is not well studied.

In Figure 5, we show the PDS of 2RXP J130159.6–635806 obtained with the NuSTAR data in the fourth observation after subtracting the pulse profile folded with the measured period from the light curve. The solid line represents the fitting model in the form of a broken power law. The measured break frequency is 0.0066 Hz (shown by dotted line), and it is clear that it is shifted toward higher frequencies relative to the spin frequency in this observation (0.0015 Hz; shown by dashed vertical line). The power-law slope above the break frequency is fixed at –2 (Revnivtsev et al. 2009). The best-fit value of the slope below the break is $-0.33$. Given the steady persistent luminosity, we conclude that the spin-up observed during the last ∼20 years is caused by a long-term mass accretion rate that is high enough to squeeze the magnetosphere inside the corotational radius. This finding confirms the uniquness of 2RXP J130159.6–635806 among the other X-ray pulsars in binary systems with Be companions.

Zoom In Zoom Out Reset image size

According to measurements from the ASCA and XMM-Newton observatories, the spectrum of 2RXP J130159.6–635806 is characterized by a moderate absorption value ${N}_{{\rm{H}}}\;=(2.48\pm 0.07)\times {10}^{22}\;{\mathrm{cm}}^{-2}$, which is stable over a dozen years (Chernyakova et al. 2005). This value was obtained by approximating the source spectra as a power-law model modified by interstellar absorption (wabs model in the XSPEC package). Note that it is just slightly higher than the value of the interstellar hydrogen absorption (1.7–1.9) × 10²² cm⁻² determined by Dickey & Lockman (1990) in the direction of 2RXP J130159.6–635806.

In the first three NuSTAR observations, 2RXP J130159.6–635806 serendipitously appeared highly offset from the optical axis, and at the large off-axis angles the absolute flux measurements can be affected by inaccurate PSF positioning and subsequent inapropriate weighting in the spectrum extraction procedure. This leads to large systematics for the extracted spectra. Therefore, we restrict detailed analysis of the source spectrum to the data obtained in the fourth (on-axis) observation, which provides high-quality data.

In general, the spectrum of 2RXP J130159.6–635806 has a shape that is typical for accreting pulsars in binary systems, showing a hard spectrum with an exponential cutoff at high energies. Figure 6(a) presents the phase-averaged spectrum approximated with the most suitable spectral model determined below. We initially modeled data with a cutoff power-law model (cutoffpl in the XSPEC package)

$$\begin{eqnarray}&&{{AE}}^{-{\rm{\Gamma }}}{e}^{\;-\;\frac{E}{{E}_{\mathrm{fold}}}},\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Gamma }}$ is the photon index, ${E}_{\mathrm{fold}}$ is the characteristic energy of the cutoff (e.g., the folding energy), and $A$ is a normalization. This model was modified by interstellar absorption in the form of the wabs model. As the working energy range of the NuSTAR observatory begins at 3 keV, it is not very sensitive to measuring low absorption columns. Therefore, in the following analysis, the interstellar absorption was fixed to the value ${N}_{{\rm{H}}}=2.48\times {10}^{22}\;{\mathrm{cm}}^{-2}$ measured by Chernyakova et al. (2005). Note that we simultaneously fitted spectra obtained by both NuSTAR modules. To take into account the uncertainty in their relative calibrations, which may be even more of a concern for the observations where the source is highly offset from the optical axis, we added a cross-calibration constant between the modules. The fitting parameters for different data sets are shown in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Parameters for the 2RXP J130159.6–635806 Phase-averaged Spectral Analysis Based on NuSTAR Observations

Seq.	Modela	Constb	Photon	${E}_{\mathrm{cut}}$	${E}_{\mathrm{fold}}$	Fluxc2–10 keV		${\chi }_{\nu }^{2}$ (dof)
Num.			Index	(keV)	(keV)	FPMA	FPMB
1	HI	$1.24\pm 0.02$	${1.40}_{-0.14}^{+0.07}$	${7.47}_{-1.55}^{+0.86}$	${17.43}_{-1.29}^{+1.79}$	${2.16}_{-0.38}^{+0.12}$	${2.68}_{-0.48}^{+0.13}$	1.06 (458)
2	HI	$0.61\pm 0.04$	1.37 (frozen)	${5.72}_{-1.90}^{+0.65}$	${15.49}_{-1.73}^{+2.62}$	${3.29}_{-1.47}^{+0.26}$	${2.01}_{-0.95}^{+0.04}$	0.92 (124)
3	HI	$0.92\pm 0.03$	1.37 (frozen)	${5.68}_{-0.62}^{+0.53}$	${17.48}_{-1.35}^{+1.45}$	${2.11}_{-0.13}^{+0.08}$	${1.95}_{-0.10}^{+0.06}$	0.81 (264)
4	CU	$1.00\pm 0.01$	$1.04\pm 0.03$	⋯	$12.62\pm 0.48$	${3.79}_{-0.10}^{+0.05}$	${3.79}_{-0.11}^{+0.05}$	1.07 (921)
	HI	$1.00\pm 0.01$	$1.32\pm 0.03$	${6.48}_{-0.15}^{+0.22}$	${16.78}_{-0.56}^{+0.68}$	${3.79}_{-0.10}^{+0.05}$	${3.79}_{-0.11}^{+0.05}$	1.01 (920)
	HI+Fed	$1.00\pm 0.01$	$1.37\pm 0.04$	${6.94}_{-0.42}^{+0.50}$	$17.96\pm 0.87$	${3.79}_{-0.10}^{+0.05}$	${3.79}_{-0.11}^{+0.05}$	1.00 (919)

Notes.

^aThe xspec spectral model used. "CU": wabs*cutoffpl, "HI": wabs*powerlaw*highecut, "HI+Fe": wabs(powerlaw*highecut+gau). ^bConstant factor of the FPMB spectrum relative to FPMA. ^cThe absorbed flux in units of ${10}^{-11}$ erg s⁻¹ cm⁻². ^dThe corresponding 6.4 keV iron line parameters are: $\sigma =0.1$ keV (fixed), normalization $(2.31\pm 0.70)\times {10}^{-5}\;\mathrm{ph}\ {\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$, and line equivalent width $\mathrm{EW}={44}_{-32}^{+48}$ eV (3σ error).

Download table as:  ASCIITypeset image

The "cutoff" model approximates the source spectrum relatively well with ${\chi }^{2}=992.18$ for 921 degrees of freedom (dof; see Table 3). Nevertheless, a wave-like structure is clearly seen in the residual panel Figure 6(b). To improve the quality of the fit, we applied another continuum model in the form of a power-law multiplied by a high-energy cutoff (highecut model in the XSPEC package; White et al. 1983). This model can be written as

$$\begin{eqnarray}{{AE}}^{-{\rm{\Gamma }}}\times \{\begin{array}{ll}1 & (E\leqslant {E}_{\mathrm{cut}})\\ {e}^{-(E-{E}_{\mathrm{cut}})/{E}_{\mathrm{fold}}} & (E\gt {E}_{\mathrm{cut}}),\end{array}\end{eqnarray} \tag{ 3 }$$

where ${E}_{\mathrm{cut}}$ is the energy where the cutoff starts. As in the previous case, we fixed the interstellar absorption value. Residuals of modeling the source spectrum with "highecut" are presented in Figure 6(c), and best-fit parameters are listed in Table 3. The "highecut" model significantly improves the quality of the fit (${\chi }^{2}=929.72$ for 920 dof). The high-energy cutoff value ${E}_{\mathrm{cut}}$ is found to be ${6.48}_{-0.15}^{+0.22}$ keV, which is significantly lower than the value of ∼25 keV reported by Chernyakova et al. (2005) using simultaneous XMM-Newton and INTEGRAL data. The apparent discrepancy is probably due to the lower statistical quality of the INTEGRAL data and the gap between the energy bands covered by the XMM-Newton and INTEGRAL observatories. An additional possible problem is that the cutoff energy ${E}_{\mathrm{cut}}=6.48$ keV is very close to the energy of the iron fluorescent line at 6.4 keV. The simplistic "highecut" model might hide the presence of the emission line in the source spectrum. Observed deviations of the measured spectrum from the model near this energy argue in favor of this possibility.

To investigate this issue, we added a Gaussian emission line to the "highecut" model, fixing its energy to 6.4 keV and width to 0.1 keV, allowing its normalization to be a free parameter. This resulted in an additional improvement of the fit to ${\chi }^{2}=916.97$ for 919 dof for a normalization of $(2.31\pm 0.70)\times {10}^{-5}\;\mathrm{ph}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. The corresponding equivalent width of the line is $\mathrm{EW}={44}_{-32}^{+48}$ eV (3σ error). We determined the significance of the line using the XSPEC script simftest with $4\;\times \;{10}^{4}$ trials, and found that presumption against the null hypothesis, or no line required by the data is $3\;\times \;{10}^{-3}$, which corresponds to $\sim 3\sigma$ line detection, assuming a normal distribution. The residuals of the "highecut" model with the 6.4 keV iron line are shown in Figure 6(d). The model itself, together with spectral data points, is shown in Figure 6(a).

Considering the moderate spectral resolution of the NuSTAR observatory and the low significance of the detected iron line, we investigated the possibility of the presence of an iron fluorescence line in the XMM-Newton spectrum of 2RXP J130159.6–635806. As mentioned above, the X-ray pulsar 2RXP J130159.6–635806 was observed with XMM-Newton many times during programs studying PSR B1259–63. For our purposes, we chose the two observations with the longest exposures—ObsID 0092820301 (∼41.2 ksec) and ObsID 0504550601 (∼55.3 ksec). As for the NuSTAR observations, we modeled the spectra of 2RXP J130159.6–635806, including the Gaussian line at 6.4 keV. For both of the XMM-Newton observations, we did not find a significant improvement in the fits when the iron line was added and obtained a conservative upper limit (3σ) for the equivalent width of the iron line of 110 eV, which is consistent with the NuSTAR results.

Despite the fact that the first three NuSTAR observations were made at large offset angles, and the corresponding statistics are significantly lower than for the fourth observation, some useful spectral information can still be extracted. For these observations, we used the same "highecut" model. Due to low statistics and poor fit, we fixed power-law slopes in the second and third data sets at a ${\rm{\Gamma }}=1.37$ value measured in the fourth observation. As seen from Table 3, where the best-fit parameters are listed, the principal parameters—power-law slope (${\rm{\Gamma }}$), cut-off energy (${E}_{\mathrm{cut}}$), and folding energy (${E}_{\mathrm{fold}}$)—of the first three high-offset observations are in general agreement with the fourth (on-axis) observation modeled with "highecut."

The estimated 2–10 keV flux of 2RXP J130159.6–635806 is about $3\;\times \;{10}^{-11}$ erg s⁻¹ cm⁻² during our observations, which is in agreement with the value of (2–3) $\times \;{10}^{-11}$ erg s⁻¹ cm⁻² measured by Chernyakova et al. (2005). An absence of strong transient activity from 2RXP J130159.6–635806 on long timescales is confirmed with the RXTE/ASM and Swift/BAT instruments. Note that Chernyakova et al. (2005) reported on a flaring episode, but the flux of the flare was relatively low, rising by a factor of a few (up to $\sim {10}^{-10}$ erg s⁻¹ cm⁻²).

Finally, in order to check for the possible presence of a cyclotron absorption line in the source spectrum, we modified the best-fit model (wabs*(powerlaw*highecut+gau)) by including an absorption feature in the form of a Lorentzian optical depth profile (cyclabs model in XSPEC; Mihara et al. 1990). The search procedure was performed following the prescription from Tsygankov & Lutovinov (2005). Namely, we varied the energy of the line over the range between 5 and 50 keV with 3 keV steps and the line width between 2 and 12 keV with 2 keV steps, leaving the line depth as a free parameter. We did not find strong evidence for a cyclotron line in the spectrum of 2RXP J130159.6–635806 (no trials gave a significance higher than $\sim 2\sigma$).

In order to study the evolution of the source spectrum over the pulse period, we performed pulse phase-resolved spectroscopy using the data from the fourth observation. The period was divided into 10 phase bins with the zero phase coinciding with the main minimum in the pulse profile. According to the pulse phase-averaged spectral analysis, the fitting model was chosen in the form of (wabs*(powerlaw*highecut+gau)). However, due to much lower statistics, the photon index was fixed at the value from the averaged spectrum (${\rm{\Gamma }}=1.37$). We selected this parameter due to its virtual constancy over the pulse in our preliminary analysis of the same data. The ${N}_{{\rm{H}}}$ value is not very well constrained by the NuSTAR data; however we checked that it is consistent with being constant in phase, and we fixed it as well.

The results are shown in Figure 7. The average pulse profile of 2RXP J130159.6–635806 across the entire NuSTAR energy range is presented in the upper panel. The lower panels demonstrate the behavior of the two free spectral parameters: the cut-off energy (${E}_{\mathrm{cut}}$) and the folding energy (${E}_{\mathrm{fold}}$). The cut-off energy is quite stable over the pulse except at the pulse minimum where its value increases approximately by a factor of two. The folding energy demonstrates an apparent correlation with the pulse intensity, which is probably caused by increasing a spectral hardness around the pulse maximum. Such behavior of the spectral parameters over the pulse can explain the corresponding behavior of the hardness ratios constructed from the pulse profiles in different energy bands, previously shown in Figure 3. The (20–40)/(10–20) keV hardness ratio demonstrates correlation with the pulse intensity and the (10–20)/(3–10) keV ratio peaks at the pulse minimum right in place where ${E}_{\mathrm{cut}}$ shifts from ∼7 to ∼16 keV. Finally, it is necessary to note that the observed behavior of spectral parameters with the pulse phase can be caused by both physical and artificial reasons (in particular, due to limitations of available data, the adopted spectral model, and which model parameters were fixed). More data are required to confidently constrain all of the spectral parameters and trace their behavior with the pulse phase.

Zoom In Zoom Out Reset image size