Spectral and Timing Analysis of the Accretion-powered Pulsar 4U 1626–67 Observed with Suzaku and NuSTAR

NuSTAR is NASA's 11th small explorer (SMEX) mission (Harrison et al. 2013). The satellite covers an energy range of 3–79 keV with two CdZnTe detectors, Focal Plane Module A and B (FPMA and FPMB), located at the focal planes of its hard X-ray mirrors. NuSTAR observed 4U 1626−67 from 2015 May 4, 12:46 UT, to May 5, 20:27 UT (MJD 57146.5319 to 57147.8521). We reduced the data with version 1.6.0 of the nupipeline software as distributed with HEASOFT 6.19. After standard screening of the data with v20170503 of the NuSTAR calibration data, the net exposure time of the observation was 65 ks. A barycentric correction was applied to the arrival time of each event in the FPMA and FPMB with the barycorr tool of HEASOFT. We extracted source spectra from FPMA and FPMB using a circle of radius 100'' centered on the source.

For the timing analysis, the background light curves were extracted from a circular region with a radius of 100'' at the corner of the NuSTAR field of view opposite to the source. For the spectral analysis, we modeled background spectra applying the NuSTAR background-fitting and -modeling tool, nuskybgd (Wik et al. 2014),¹⁵ using three blank-sky spectra for each telescope extracted from annular regions with radii 300''–400'', 400''–500'', and 500''–600''.

Suzaku was the fifth Japanese X-ray satellite (Mitsuda et al. 2007). It was equipped with two types of instruments, the X-ray Imaging Spectrometer (XIS; Koyama et al. 2007) and the Hard X-ray Detector (HXD; Takahashi et al. 2007). The XIS was a set of four charge-coupled device cameras at the foci of four X-ray telescopes (Serlemitsos et al. 2007), covering the lower energy band of 0.2–12 keV. XIS 0, XIS 2, and XIS 3 were front-illuminated devices, while XIS 1 was back illuminated. XIS 2 was taken offline after micrometeorite strikes prior to these observations, so we use XIS 0 and XIS 3 (combined: XIS-FI) as well as XIS 1 (XIS-BI) in this analysis. The HXD consisted of PIN silicon diodes (HXD-PIN) and Gd₂SiO₅Ce (GSO) crystal scintillators and covered the 10–600 keV energy band.

The data used here are from the second Suzaku observation of 4U 1626−67, which was performed between 2010 September 6, 12:59 UT, and 2010 September 7, 05:42 UT (MJD 55445.5410 to 55446.2375). The first observation of the source had been conducted in 2006 March and was analyzed by Iwakiri et al. (2012). The XIS was operated with the 1/4 window option and standard clocking, which has a time resolution of 2 s. As this is comparable to the 7.7 s pulse period of 4U 1626−67, phase-resolved analysis of the Suzaku data is not practical. Therefore, we used the Suzaku data for the phase-averaged analysis only. We reprocessed and screened the XIS and HXD data with the standard criteria, using the Suzaku reprocessing tool aepipeline in the HEASOFT v6.19 package with calibration versions hxd-20110913, xis-20151005, and xrt-20110630. The XIS spectra and light curves were extracted from a circular region with a radius of 43 centered at the source. We accumulated the XIS background spectra from a source-free region. We checked the pile-up effect in the same way as in Yamada et al. (2012) and found that the pile-up fraction was 3% and 1% at 5.7 and 28.7 pixels from the center of the image, respectively. We excluded the regions that show >1% pile-up fraction so that the pile-up effect is negligible for the subsequent analysis. In this paper, we only used the XIS data down to 3 keV to adjust the energy range of NuSTAR.

Because the HXD was not capable of imaging, we applied the simulated non-X-ray background model provided by the Suzaku HXD team (Fukazawa et al. 2009) and a cosmic X-ray background, assumed to be the same as the typical model obtained by HEAO-1 (Boldt 1987). The net exposures were 20.1 ks for the XIS 0, 1, 3 and 18.7 ks for the HXD, respectively. Different screening criteria for South Atlantic Anomaly passages and Earth elevation angles result in different exposures for the XIS and the HXD.

Figure 1 shows the long-term light curves and hardness ratio of 4U 1626−67 from 2009 October to 2016 January, obtained with MAXI (Matsuoka et al. 2009) on board the International Space Station (ISS). To avoid any systematic modulation with a period of 70 days caused by the ISS orbit, the data are binned to a resolution of 70 days per bin. Over the period considered here, the X-ray fluxes monotonically increased by ∼20% in both the 2–10 and 10–20 keV bands, while the hardness ratio remained constant.

Zoom In Zoom Out Reset image size

Applying epoch folding (Leahy et al. 1983) to the light curves obtained with the Suzaku HXD-PIN and NuSTAR, we find the spin period of 4U 1626−67 to be P_Suzaku = 7.6774(1) s and P_NuSTAR = 7.67295(1) s, respectively. These results are consistent with Fermi/gamma-ray burst monitory (GBM) results.¹⁶

We estimate the period derivatives, $\dot{P}$, at the epochs of the Suzaku and NuSTAR observations by performing a linear fit to the Fermi/GBM data for 60 days before and after each observation. The results are ${\dot{P}}_{{Suzaku}}=-2.8\times {10}^{-11}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ and ${\dot{P}}_{{NuSTAR}}\,=-3.3\times {10}^{-11}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, respectively. Using the conservative approach of Takagi et al. (2016), the uncertainty of these $\dot{P}$ values is 0.2 × 10⁻¹¹ s s⁻¹. Thus, the period derivative $\dot{P}$, along with the X-ray flux (Figure 1), increased significantly between the two epochs of the Suzaku and NuSTAR observations. Accretion torque theory (Ghosh & Lamb 1979) implies that this change in $\dot{P}$ is due to an increase in the accretion rate on the neutron star, consistent with the study of Takagi et al. (2016), who find from observations of 4U 1626−67 spanning 30 yr that the period derivative changes are in good agreement with the prediction by Ghosh & Lamb (1979).

3.2.1. Crab Ratio

First, we examine the ratio of the spectra to that of the Crab Nebula and pulsar in order to get a model-independent view of their overall properties. Figure 2 shows the ratio of the phase-averaged spectra obtained by Suzaku and NuSTAR to that of the Crab, using a canonical model based on the Suzaku calibration results¹⁷ as a template for the Crab spectrum. We fixed the cross-normalization factor between XIS 0 and the HXD-PIN to the standard value of 1.16. According to the NuSTAR calibration results of Madsen et al. (2017), the cross-normalization factor between the FPMA and Suzaku's instruments are 0.99 ±0.01 for XIS 0 and 1.13 ± 0.03 for the HXD-PIN, respectively. Comparing with these results, the cross-normalization uncertainty is at the few-percent level.

Zoom In Zoom Out Reset image size

The CRSF is clearly visible at ∼40 keV. The ratios in the soft X-rays below 20 keV show that the X-ray flux obtained with NuSTAR is brighter than that obtained with Suzaku, which is consistent with the long-term X-ray variation observed by MAXI (Figure 1). Above 30 keV, however, the hard X-ray flux decreased in the latter observation. In addition, the spectral peak energy appears to have slightly changed from ∼22.5 keV in the 2010 Suzaku observation to ∼21.5 keV in the 2015 NuSTAR observation. A similar spectral difference is also reported by D17, who compared the long-term Swift/Burst Alert Telescope spectrum with the NuSTAR data. D17 suggest that the difference is due to a small change of the depth of the fundamental CRSF; however, the model-independent Crab ratios in Figure 2 indicate that changes in the continuum spectral shape should also be considered.

3.2.2. Continuum

Next we modeled the phase-averaged spectra using the negative and positive power law times exponential (NPEX) model (Mihara 1995; Makishima et al. 1999). To model the CRSF, we multiplied the continuum by a Gaussian optical depth absorption model (GABS in XSPEC) $\exp (-S(E))$, where

$$\begin{eqnarray}&&S(E)=\displaystyle \frac{{d}_{{\rm{a}}}}{\sqrt{2\pi }{\sigma }_{{\rm{a}}}}{e}^{-{\left(E-{E}_{{\rm{a}}}\right)}^{2}/(2{\sigma }_{{\rm{a}}}^{2})},\end{eqnarray} \tag{ 1 }$$

and where E_a, σ_a, and d_a are the energy, width, and line depth. We also apply an iron line modeled by an additive Gaussian, and an additional blackbody for the soft excess (Schulz et al. 2001) modeled by a BBODYRAD model, to the Suzaku and NuSTAR data sets separately. The hydrogen column density is fixed at the total Galactic H i column density (${N}_{{\rm{H}}}\,=1\times {10}^{21}\,{\mathrm{cm}}^{-2}$).¹⁸

The NPEX model succeeded in reproducing the Suzaku data (left of Figure 3(a), with χ²/dof = 326.08/317) but failed to reproduce the NuSTAR data (right of Figure 3(a), with χ²/dof = 1005.1/427). However, we note that it can describe the NuSTAR spectrum quite well if the energies around the spectral peak (between 10 and 30 keV) are excluded. Figure 3(b) displays the results of a fit when the 10–30 keV data are ignored. This motivates the construction of a new model, which we term "modNPEX," consisting of an NPEX continuum with an additional broad Gaussian emission feature around 20 keV. Similar bumps have been seen in other accretion-powered pulsars (Coburn et al. 2002; Klochkov et al. 2007; Ferrigno et al. 2009; Vasco et al. 2013; Farinelli et al. 2016; Ballhausen et al. 2017). When fitting the Suzaku data, the center energy and the width of the Gaussian model were fixed at the best-fit values of the NuSTAR data. We found that both data sets are well reproduced by this model (Figure 3(c)). A broad Fe Kα emission line is also detected at ∼6.76 keV with a σ of ∼150 eV in both the observations. These line centroid energies and widths are consistent with the results in previous papers (Camero-Arranz et al. 2012, Suzaku data; D17, NuSTAR data). We also tried different phenomenological continuum models, the Fermi–Dirac cutoff model (Tanaka 1986) and a power law modified with a high-energy cut-off (HIGHECUT) model (White et al. 1983; Coburn et al. 2002; Fürst et al. 2013). As a result, we found that the modNPEX is the most successful in providing a good and consistent description of both the phase-averaged and the two phase-resolved spectra (Section 3.4). Therefore, we adopted only the modNPEX as the continuum emission model in this paper.

Zoom In Zoom Out Reset image size

It should be noted here that D17, using a HIGHECUT model (White et al. 1983; Coburn et al. 2002; Fürst et al. 2013) or bulk+thermal comptonization continuum model (Becker & Wolff 2007), noted absorption-like residuals around ∼60 keV. They therefore added a second Gaussian absorption model to describe this feature, which they interpret as the harmonic of the CRSF. However, as D17 acknowledge, the centroid energy of the feature deviates significantly from an integer multiple of the fundamental energy, quite a bit more than one expects for a harmonic CRSF (even if relativistic corrections are taken into account, the harmonic should lie at roughly twice the energy of the fundamental; see Pottschmidt et al. 2005). Additionally, the line parameters are model dependent: D17 find a significantly lower energy of 61 ± 1 keV when using the Becker & Wolff (2007) physical continuum model, compared with their HIGHECUT fits (67 ± 3 keV), while the line depth drops from ${50}_{-10}^{+17}$ to 22 ± 5.

In our analysis of the NuSTAR spectra, a second GABS feature is not needed to obtain an acceptable fit. Additionally, even when residuals around 60 keV are present, they are limited to a single bin at the extreme upper end of the useful NuSTAR spectrum. Thus, the feature reported by D17 is mainly constrained by the Swift-BAT spectrum. We further note that the Swift-BAT has a tagged ²⁴¹Am calibration source on board, which emits 59 keV photons (Gehrels et al. 2004); this contributes to its background and could result in line-like features. Thus, further study of 4U 1626−67 with energy coverage considerably higher than 60 keV is needed to properly study any possible harmonic CRSF.

As the Crab ratios from Section 3.2.1 suggest, our spectral modeling shows a decrease in the pseudo-plasma temperature kT from ${7.1}_{-0.4}^{+0.6}\,\mathrm{keV}$ to 5.9 ± 0.1 keV, while the CRSF is narrower and shallower in the NuSTAR observation compared with Suzaku. The unabsorbed X-ray fluxes in the 3.0–10.0, 30.0–60.0, and 3.0–60.0 keV bands derived with the modNPEX model fitting of the Suzaku data are ${3.21}_{-0.12}^{+0.09}\times {10}^{-10}$, ${1.72}_{-0.28}^{+0.31}\times {10}^{-10}$, and ${1.29}_{-0.14}^{+0.14}\times {10}^{-9}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, respectively. Assuming a distance of 10 kpc (optical observations constrain its distance to 5–13 kpc; Chakrabarty 1998), these fluxes correspond to an X-ray luminosity of ${1.54}_{-0.17}^{+0.17}\,\times {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the 3.0–60.0 keV band. The NuSTAR fluxes in the same energy ranges are ${3.44}_{-0.06}^{+0.04}\times {10}^{-10}$, ${1.14}_{-0.08}^{+0.07}\,\times {10}^{-10}$, and ${1.25}_{-0.04}^{+0.05}\times {10}^{-9}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, corresponding to an X-ray luminosity of ${1.50}_{-0.05}^{+0.06}\times {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the 3.0–60.0 keV band. Comparing the observed X-ray fluxes by Suzaku with the NuSTAR results, we found that the 3.0–10.0 keV flux increased by ∼5%, whereas the 30.0–60.0 keV flux decreased by ∼35%.

Table 1.  Best-fit Parameters for the Phase-averaged Suzaku and NuSTAR Spectra of 4U 1626−67 Using the modNPEX Continuum Models with One GABS Component to Describe the Fundamental Cyclotron Line

	modNPEX × GABS
	Suzaku	NuSTAR
kTBB (keV)	${0.47}_{-0.06}^{+0.06}$	${0.47}_{-0.04}^{+0.04}$
ABBa	${258}_{-104}^{+271}$	${254}_{-85}^{+174}$
α	${0.53}_{-0.18}^{+0.16}$	${0.47}_{-0.10}^{+0.09}$
kTNPEX (keV)	${7.1}_{-0.4}^{+0.6}$	${5.9}_{-0.1}^{+0.1}$
An(×10−2)b	${2.3}_{-0.6}^{+0.7}$	${2.6}_{-0.4}^{+0.4}$
Ap(×10−5)b	${4.0}_{-1.8}^{+1.4}$	${7.0}_{-1.3}^{+0.9}$
EGaussian (keV)	19.8(fix)	${19.8}_{-0.5}^{+0.4}$
σGaussian (keV)	4.8(fix)	${4.8}_{-0.5}^{+0.6}$
AGaussian(×10−3)d	${1.8}_{-1.7}^{+2.3}$	${3.1}_{-0.8}^{+1.2}$
EFe (keV)	${6.75}_{-0.06}^{+0.06}$	${6.76}_{-0.05}^{+0.05}$
σFe (keV)	${0.13}_{-0.08}^{+0.06}$	${0.15}_{-0.06}^{+0.06}$
AFe(×10−4)c	${1.3}_{-0.4}^{+0.5}$	${1.2}_{-0.2}^{+0.3}$
ECRSF (keV)	${38.2}_{-0.8}^{+0.9}$	${37.7}_{-0.1}^{+0.1}$
σCRSF (keV)	${5.1}_{-0.7}^{+0.7}$	${4.2}_{-0.1}^{+0.1}$
dCRSF	${22.9}_{-2.6}^{+5.7}$	${14.7}_{-0.8}^{+0.8}$
${C}_{{\mathrm{XIS}}_{\mathrm{BI}}}$ d	${0.941}_{-0.007}^{+0.007}$	⋯
CFPMBe	⋯	${1.002}_{-0.002}^{+0.002}$
χ2/dof	323.48/316	533.08/422
pnullf	3.7 × 10−1	1.9 × 10−4

Notes.

^a ${R}_{\mathrm{km}}^{2}/{d}_{10}^{2}$, where R_km is the radius of the blackbody in kilometers and d₁₀ is the distance in units of 10 kpc. ^bNormalization of the power law, defined at 1 keV in units of photons keV⁻¹ cm⁻² s⁻¹. ^cNormalization of the Gaussian, defined in units of photons keV⁻¹ cm⁻² s⁻¹. ^dA cross-calibration constant for XIS-BI relative to XIS-FI. ^eA cross-calibration constant for FPMB relative to FPMA. ^fNull hypothesis probability.

Download table as:  ASCIITypeset image

3.2.3. An Asymmetric Cyclotron Line

The residuals of the modNPEX fit to the NuSTAR data show a slight systematic structure remaining at around the fundamental CRSF energy (Figure 3(c), right). The presence of a complex line profile for the fundamental CRSF is consistent with D17, despite their different approach to modeling the spectrum. We tried different models to reproduce this feature better and investigate its shape. First, we replaced the GABS model by a pseudo-Lorentzian absorption model (CYAB model, cyclabs in XSPEC) and fitted the spectrum using modNPEX as a continuum model. As can be seen in Figure 4 and Table 2, this model did not improve the fit relative to the GABS CRSF model. Next, we added another GABS model, as has been used for V 0332+53 (X 0331+53) (Kreykenbohm et al. 2005; Pottschmidt et al. 2005; Nakajima et al. 2010) and Cep X-4 (Fürst et al. 2015). With this model, we successfully eliminated the residuals around the CRSF (Figure 4 and Table 2). These results are largely consistent with D17, although they tie the line energies of the GABS models to having the same values.

Zoom In Zoom Out Reset image size

Table 2.  Best-fit Parameters for the Phase-averaged NuSTAR Spectra of 4U 1626−67 Using the modNPEX Continuum Models with Different One- and Two-component Models to Describe the Fundamental Cyclotron Line

	modNPEX	modNPEX	(modNPEX + Gaussian)
	× CYAB	× GABS × GABS	× GABS
kTBB (keV)	${0.47}_{-0.04}^{+0.04}$	${0.47}_{-0.03}^{+0.03}$	${0.47}_{-0.03}^{+0.03}$
ABBa	${263}_{-89}^{+167}$	${266}_{-86}^{+154}$	${266}_{-87}^{+155}$
α	${0.52}_{-0.07}^{+0.06}$	${0.48}_{-0.05}^{+0.05}$	${0.48}_{-0.05}^{+0.05}$
kTNPEX (keV)	${6.3}_{-0.1}^{+0.1}$	${5.9}_{-0.1}^{+0.1}$	${5.9}_{-0.1}^{+0.1}$
An(×10−2)b	${2.6}_{-0.4}^{+0.3}$	${2.6}_{-0.2}^{+0.2}$	${2.6}_{-0.2}^{+0.2}$
Ap(×10−5)b	${6.1}_{-1.0}^{+0.7}$	${7.1}_{-0.4}^{+0.4}$	${7.1}_{-0.4}^{+0.4}$
EGaussian (keV)	${20.6}_{-0.6}^{+0.4}$	19.8(fix)	19.8(fix)
${\sigma }_{\mathrm{Gaussian}}$ (keV)	${5.0}_{-0.4}^{+0.4}$	4.8(fix)	4.8(fix)
AGaussian(×10−3)c	${3.1}_{-0.8}^{+1.3}$	${3.1}_{-0.1}^{+0.3}$	${3.1}_{-0.1}^{+0.3}$
EFe (keV)	${6.76}_{-0.05}^{+0.05}$	${6.76}_{-0.05}^{+0.05}$	${6.76}_{-0.05}^{+0.05}$
σFe (keV)	${0.15}_{-0.06}^{+0.06}$	${0.15}_{-0.06}^{+0.06}$	${0.15}_{-0.06}^{+0.06}$
AFe(×10−4)c	${1.2}_{-0.2}^{+0.3}$	${1.2}_{-0.2}^{+0.3}$	${1.2}_{-0.2}^{+0.3}$
ECRSF (keV)	${36.6}_{-0.1}^{+0.1}$	${39.9}_{-1.3}^{+0.6}$	${39.2}_{-0.7}^{+0.8}$
σCRSF (keV)	${5.3}_{-0.3}^{+0.3}$	${2.6}_{-0.8}^{+0.6}$	${2.9}_{-0.5}^{+0.5}$
τCRSF	${1.6}_{-0.1}^{+0.1}$	⋯	⋯
dCRSF	⋯	${5.3}_{-3.8}^{+5.9}$	${8.8}_{-3.2}^{+2.8}$
Eabs (keV)	⋯	${35.9}_{-2.7}^{+1.4}$	${32.8}_{-1.1}^{+0.9}$
σabs (keV)	⋯	${3.6}_{-0.7}^{+0.5}$	${3.4}_{-0.5}^{+0.5}$
dabs	⋯	${8.7}_{-5.6}^{+4.2}$	⋯
Aabs(×10−3)c	⋯	⋯	$-{1.2}_{-0.4}^{+0.4}$
CFPMBd	${1.002}_{-0.002}^{+0.002}$	${1.002}_{-0.002}^{+0.002}$	${1.002}_{-0.002}^{+0.002}$
χ2/dof	596.91/422	508.08/421	507.98/421
pnulle	2.1 × 10−9	2.3 × 10−3	2.3 × 10−3

Notes.

^a ${R}_{\mathrm{km}}^{2}/{d}_{10}^{2}$, where R_km is the radius of the blackbody in kilometers and d₁₀ is the distance in units of 10 kpc. ^bNormalization of the power law, defined at 1 keV in units of photons keV⁻¹ cm⁻² s⁻¹. ^cNormalization of the Gaussian, defined in units of photons keV⁻¹ cm⁻² s⁻¹. ^dA cross-calibration constant for FPMB relative to FPMA. ^eNull hypothesis probability.

Download table as:  ASCIITypeset image

However, because the GABS model only produces an absorption feature, we also tried replacing the second GABS component with an additive Gaussian model (i.e., [modNPEX +Gaussian] × GABS), allowing the Gaussian normalization to be both positive and negative. The best-fit normalization using the model is negative at the 5σ level, and the energy was consistent with that found using the 2 × GABS CRSF. To evaluate the chance probability of improvement of adding the extra Gaussian component, we simulated 400,000 data sets using simftest in XSPEC. The evaluated chance probability was 6.3 × 10⁻⁵. Therefore, we have possibly detected a complex line profile for the fundamental CRSF in the phase-averaged spectrum.

Because the energy spectrum of 4U 1626−67 is known to depend strongly on pulse phase (Pravdo et al. 1979; Kii et al. 1986; Iwakiri et al. 2012), we present in these next two sections a detailed phase-resolved study. First, to investigate the geometry of the neutron star quantitatively, we performed pulse profile modeling using a new relativistic ray tracing code. Due to the low time resolution of the Suzaku data below 10 keV, we concentrate on the NuSTAR data. Figure 5 shows the energy-resolved and background-subtracted pulse profiles in 11 different energy bands obtained with NuSTAR. The pulse profiles strongly depend on energy, with a double-peaked structure below 10 keV, which becomes single peaked and almost sinusoidal in the higher energy bands. The characteristics of the pulse profiles seen by NuSTAR are consistent with the RXTE results observed in 2010 (Beri et al. 2014).

Zoom In Zoom Out Reset image size

We use the relativistic light-bending code of Falkner et al. (Falkner et al. 2018a, 2018b, see also Schönherr et al. 2014) to model the energy-resolved pulse profiles of 4U 1626−67 (see Figure 5). This code follows a similar approach to that of Ferrigno et al. (2011) to obtain the observable energy and phase-dependent flux. In contrast to Ferrigno et al. (2011), we are able to apply arbitrary emission patterns to emission regions of any geometrical shape. The code has been used previously by Schönherr et al. (2014) to investigate the energy-dependent phase lags in accreting neutron stars.

In our model for 4U 1626−67 we consider a canonical neutron star of mass M = 1.4 M_⊙ and radius R = 10 km. The observed X-rays are emitted by two cylindrical accretion columns AC1 and AC2 of height ${h}_{\mathrm{AC}1,\mathrm{AC}2}$ and radius ${r}_{\mathrm{AC}1,\mathrm{AC}2}$. Allowing for an asymmetric magnetic field, the columns are positioned individually at azimuthal angles ${{\rm{\Phi }}}_{\mathrm{AC}1,\mathrm{AC}2}$ and polar angles ${{\rm{\Theta }}}_{\mathrm{AC}1,\mathrm{AC}2}$, respectively, in a coordinate system that is measured with respect to the neutron star's rotational axis. The angle between the line of sight and the neutron star's angular momentum vector specifies the inclination i of the neutron star. Hence i = 0° would correspond to a face-on system.

We make the simplified assumption that the emission pattern of the columns can be described as a mixture of Gaussian-like fan- and pencil-beam emission components in the frame of rest of the neutron star's surface, e.g., at a given energy the emissivity of one accretion column is given by

$$\begin{eqnarray}&&{I}_{E}(\eta )={N}_{{\rm{p}}}\exp \left(-{\left[\displaystyle \frac{\eta -{\bar{\eta }}_{{\rm{p}}}}{\sqrt{2}{\sigma }_{{\rm{p}}}}\right]}^{2}\right)+{N}_{{\rm{f}}}\exp \left(-{\left[\displaystyle \frac{\eta -{\bar{\eta }}_{{\rm{f}}}}{\sqrt{2}{\sigma }_{{\rm{f}}}}\right]}^{2}\right),\end{eqnarray} \tag{ 2 }$$

where η is the angle of the emitted photons measured with respect to the magnetic field axis in the frame of rest of the neutron star surface, and where the energy dependent quantities $\bar{\eta }$, σ, and N describe the direction of peak emissivity, the width, and the strength of the pencil-beam (p) and fan-beam (f) components, respectively. For the fan beam we set ${\bar{\eta }}_{{\rm{f}}}=90^\circ$, i.e., the fan beam is fixed to emit from the sides of the accretion column perpendicular to the B-field. While we impose the same model for the emissivity pattern on both poles, the fluxes of the beams are allowed to vary freely.

We model the energy-dependence of the pulse profile by letting the parameters N, σ, and $\bar{\eta }$ of the emission pattern be energy dependent. For the models described in the following we assume that the emissivity of the accretion column is independent of height and thus constant over the whole column. From Equation (2) we then derive the observed energy- and phase-dependent total flux

$$\begin{eqnarray}&&{{ \mathcal F }}_{E}(\phi )={F}_{1}(\phi ,{I}_{E})+{F}_{2}(\phi ,{I}_{E}),\end{eqnarray} \tag{ 3 }$$

where ϕ is the pulse phase and F_1,2 is the flux of the individual accretion columns emitting with the given, energy-dependent, emission pattern I_E (Equation (2)). The description of the calculation of F_1,2, which accounts for all general relativistic effects and is in addition to the parameters here dependent on the neutron star's inclination, i, is beyond the scope of this paper and is given in Falkner et al. (2018a, 2018b).

Figure 5 includes the best-fit model described above for the energy-resolved and background-subtracted NuSTAR pulse profiles in 11 energy bands, which were fitted simultaneously. Lines show the overall model and the individual contributions of the fan and pencil beams from each accretion column. The corresponding parameters are listed in Table 3. The geometrical parameters (i, ${{\rm{\Theta }}}_{\mathrm{AC}1,\mathrm{AC}2}$, ${{\rm{\Phi }}}_{\mathrm{AC}1,\mathrm{AC}2}$, ${h}_{\mathrm{AC}1,\mathrm{AC}2}$, ${r}_{\mathrm{AC}1,\mathrm{AC}2}$) are global, that is, the same for all energies, while the parameters describing the emission profile are determined for each individual energy band. In the best-fit case the size of the second accretion column is tied to the first one, i.e., ${r}_{\mathrm{AC}2}={r}_{\mathrm{AC}1}$ and ${h}_{\mathrm{AC}2}={h}_{\mathrm{AC}1}$. Further, the shape of the energy-dependent emission profiles of the two columns are tied together, i.e., ${\sigma }_{{{\rm{p}}}_{2}}={\sigma }_{{{\rm{p}}}_{1}}$, ${\bar{\eta }}_{{{\rm{p}}}_{2}}={\bar{\eta }}_{{{\rm{p}}}_{1}}$ for the pencil-beam component and ${\sigma }_{{{\rm{f}}}_{2}}={\sigma }_{{{\rm{f}}}_{1}}$, ${\bar{\eta }}_{{{\rm{f}}}_{2}}={\bar{\eta }}_{{{\rm{f}}}_{1}}=90^\circ$ for the fan-beam component. We found that there are two local minima of the χ² landscape of the observer inclination i; one solution is i = 9°, and the other is i = 27°. Because the 9° solution shows a simpler pulse profile evolution than the i = 27° solution and because of its consistency with the physical simulation results taking into account the anisotropy of the scattering cross section in a strong magnetic field according to Kii et al. (1986) (see details in Section 4.3), we only show the results from the i = 9° solution in this paper. Table 3 also gives the 90% confidence levels for the parameters. These uncertainties are purely statistical and driven by a complex χ² landscape. It is probable that these uncertainties are systematically underestimated. In relation to each other the uncertainties, however, indicate that the radius and height of the column, and for some energies the width of the fan beam, are much less constrained than the other parameters.

Table 3.  Fit Parameters of the Accretion Column Model

Global	Column1					Column2
i	ΦAC1	ΘAC1	rAC1	${h}_{\mathrm{AC}1}$		ΦAC2	ΘAC2	χ2/dof
(deg)	(deg)	(deg)	(m)	(m)		(deg)	(deg)
9.5 ± 0.2	73.6 ± 0.3	${12.50}_{-0.11}^{+0.07}$	590 ± 50	${252}_{-14}^{+12}$		290 ± 2	158.1 ± 0.2	1146/631
	pencil1			fan1		fan2
ΔE	${N}_{{{\rm{p}}}_{1}}$	${\sigma }_{{{\rm{p}}}_{1}}$	${\bar{\eta }}_{{{\rm{p}}}_{1}}$	${N}_{{{\rm{f}}}_{1}}$	${\sigma }_{{{\rm{f}}}_{1}}$	${N}_{{{\rm{f}}}_{2}}$
(keV)	(cts s−1)	(deg)	(deg)	(cts s−1)	(deg)	(cts s−1)
3–5	${0.95}_{-0.02}^{+0.01}$	${4.21}_{-0.06}^{+0.07}$	5.73 ± 0.06	${1.21}_{-0.02}^{+0.03}$	${180}_{-16}^{+0}$	${0.011}_{-0.011}^{+0.013}$
5–7	1.14 ± 0.02	4.37 ± 0.06	6.12 ± 0.06	${1.09}_{-0.02}^{+0.03}$	${180}_{-26}^{+0}$	0.08 ± 0.02
7–9	${1.03}_{-0.02}^{+0.03}$	4.66 ± 0.05	${6.65}_{-0.06}^{+0.07}$	${0.84}_{-0.04}^{+0.03}$	${180}_{-67}^{+0}$	0.16 ± 0.02
9–11	${0.83}_{-0.03}^{+0.06}$	${4.84}_{-0.06}^{+0.05}$	${7.33}_{-0.08}^{+0.10}$	${0.60}_{-0.06}^{+0.08}$	${61}_{-10}^{+120}$	0.17 ± 0.02
11–13	${0.77}_{-0.09}^{+0.04}$	7.3 ± 0.7	9.4 ± 0.2	${0.17}_{-0.04}^{+0.08}$	${20}_{-6}^{+12}$	0.13 ± 0.02
13–15	${0.56}_{-0.03}^{+0.02}$	7.8 ± 0.6	${10.7}_{-0.2}^{+0.3}$	${0.11}_{-0.02}^{+0.03}$	${14.8}_{-0.7}^{+4.0}$	0.11 ± 0.02
15–17	${0.46}_{-0.02}^{+0.03}$	8.4 ± 0.5	12.4 ± 0.4	0.06 ± 0.02	${13.0}_{-1.7}^{+0.9}$	0.082 ± 0.009
17–20	${0.55}_{-0.04}^{+0.02}$	8.3 ± 0.3	${14.8}_{-0.6}^{+0.3}$	${0.026}_{-0.008}^{+0.050}$	${7}_{-3}^{+6}$	0.08 ± 0.01
20–25	0.53 ± 0.01	7.8 ± 0.2	16.2 ± 0.3	${0.024}_{-0.005}^{+0.006}$	${5.0}_{-0.3}^{+1.7}$	${0.067}_{-0.008}^{+0.007}$
25–30	${0.18}_{-0.03}^{+0.02}$	6.9 ± 0.3	${17.0}_{-0.7}^{+0.5}$	${0.02}_{-0.01}^{+0.02}$	${11}_{-5}^{+2}$	0.025 ± 0.005
30–50	${0.070}_{-0.008}^{+0.004}$	${8.3}_{-0.5}^{+0.6}$	18.1 ± 0.7	$\leqslant 0.007$	${14}_{-8}^{+167}$	0.018 ± 0.004

Note. Parameters not listed here are fixed or tied. In particular we impose r_AC1 = r_AC2 and h_AC1 = h_AC2 for the column size, ${N}_{{{\rm{p}}}_{2}}={N}_{{{\rm{p}}}_{1}}$, ${\sigma }_{{{\rm{p}}}_{2}}={\sigma }_{{{\rm{p}}}_{1}}$, ${\bar{\eta }}_{{{\rm{p}}}_{2}}={\bar{\eta }}_{{{\rm{p}}}_{1}}$ for the pencil beam, and ${\sigma }_{{{\rm{f}}}_{2}}={\sigma }_{{{\rm{f}}}_{1}}$ and ${\bar{\eta }}_{{{\rm{f}}}_{2}}={\bar{\eta }}_{{{\rm{f}}}_{1}}=90^\circ$ for the fan beam. The given errors correspond to the 90% confidence level.

Download table as:  ASCIITypeset image

The derived best-fit values of ${{\rm{\Theta }}}_{\mathrm{AC}1,\mathrm{AC}2}$, ${{\rm{\Phi }}}_{\mathrm{AC}1,\mathrm{AC}2}$ indicate an asymmetric B-field configuration of 4U 1626−67. Such an orientation of a magnetic field axis was also suggested in the previous pulse profile modeling results of Leahy (1991), which show that nine of 20 pulsars require a magnetic field axis offset. The fit shows that the magnetic field of the first column passes through close to the line of sight during each rotation. That is, at pulse phase ϕ = Φ_AC1, when the first column is in the front, we are looking at the first column from above with an angle to its magnetic field axis of approximately Θ_AC1 − i = 3°. Furthermore, the displacement of the two accretion columns compared with the symmetric antipodal case is given as

$$\begin{eqnarray}&&{\rm{\Delta }}=\pi -\arccos \left({{\boldsymbol{n}}}_{\mathrm{AC}1}\,\circ \,{{\boldsymbol{n}}}_{\mathrm{AC}2}\right),\end{eqnarray} \tag{ 4 }$$

which is the angular distance of the two unit vectors, ${{\boldsymbol{n}}}_{\mathrm{AC}1,\mathrm{AC}2}$, of the accretion columns positions corresponding to their azimuthal and polar angles. The geometry of our best-fit model yields a displacement of Δ = 139, which represents a moderate asymmetry. In our best-fit model we tie the accretion columns dimensions to prevent parameter degeneracy. Such an asymmetric B-field configuration, however, may suggest columns of different sizes due to possible asymmetric accretion flows, that is, different accretion rates for the two poles (see, e.g., Becker & Wolff 2007; Postnov et al. 2015). Therefore, and due to the fact that we make the very simple assumption of a homogeneously emitting surface, the obtained heights and radii of the columns should not be interpreted as physical quantities.

Despite the simplified assumptions entering the beam pattern, our model describes the observed pulse profiles and their energy evolution remarkably well and with a smooth variation of all relevant parameters of the emission characteristics. The evolution of the pulse profile is characterized by a very wide fan beam, which strongly decreases in width as the energy increases. In contrast, the shape of the pencil-beam component only slightly changes, but the direction of its peak emission changes significantly. This shift explains the widening of the gap between the double peaks in the pulse profile. From 3 to 10 keV the parameters corresponding to the first column change only slightly, but the normalization of the fan beam of the second column increases, explaining the changes seen in the pulse profile at these lower energies. For the geometry in our best-fit model the pencil-beam component of the second accretion column is directed away from the observer over the whole rotational phase and is therefore never observed. To ensure the correctness¹⁹ of our fit we therefore tie the normalization of the two pencil-beam components, i.e., ${N}_{{{\rm{p}}}_{2}}={N}_{{{\rm{p}}}_{1}}$. Figure 6 shows the energy dependency of the parameter values.

Zoom In Zoom Out Reset image size

The decomposition of the pulse profiles in the middle columns of Figure 5 illustrates how these parameter changes manage to reproduce the pulse profiles so well: the pencil beam is responsible for the distinct and symmetric double peak that characterizes the softer energy bands. The peaks are close together because the pencil beam is directed upward, with only a small offset, ${\bar{\eta }}_{{\rm{p}}}$, to the magnetic field, which resembles a conical emission pattern. The strong nonpulsed continuum between the two peaks is produced by the broad fan beam, which is shifted by half a phase with respect to the dip between the double peak of the pencil beam. With increasing energy, the double peak decreases in importance and its width broadens while the gap in between increases. This behavior is reflected in the best-fit parameters by showing that the offset angle ${\bar{\eta }}_{{\rm{p}}}$ and beamwidth σ_p increase with energy (see Figure 6). In addition, the flat plateau at low energies evolves into an asymmetric peak, which is caused by the narrowing of the fan beam. The asymmetry visible in this pulse is caused by the fan beam of the slightly misaligned second accretion column. A consequence of this best-fit geometry is that the pencil beam of the second column is directed away from the observer at all pulse phases, and thus it is not observable. A second consequence is that the evolution of the direction of the peak emissivity of the primary pencil beam, ${\bar{\eta }}_{{{\rm{p}}}_{1}}$, dominates the evolution of the pulsed fraction, $({{ \mathcal F }}_{\max }-{{ \mathcal F }}_{\min })/{{ \mathcal F }}_{\max }$, while the width of the fan beam, σ_f, has only a minor influence (Figure 6).

We note that extrapolating the behavior of the pencil and fan beams to even lower energies than considered here predicts an evolution of the pulse profile toward a shape dominated by the single broad hump of the fan beam, consistent with the pulse profiles seen by XMM-Newton (Krauss et al. 2007) and Chandra (Hemphill et al. 2017).

Modeling the pulse profiles also yields our viewing angle, η, onto the two accretion columns. This parameter is important for the interpretation of the CRSF, the shape of which strongly depends on the angle under which we see the magnetic field (Schwarm et al. 2017a, 2017b). As shown in Figure 7, η is strongly influenced by the effect of light bending, for the first column η varies between 5° and 22° in a small band with mean width $\sim 1^\circ$, whereas for the second column the mean width of the band is ∼11° between 91° and 110°.

Zoom In Zoom Out Reset image size

Compared with other models put forward for explaining the energy-dependent change of the pulse profile of 4U 1626−67, our pulse decomposition explains the observed energy-dependent behavior solely by a change in the emission characteristics of the accretion column, without invoking foreground effects, such as the absorption by an accretion stream proposed by Beri et al. (2014). The simpler explanation is possible by virtue of the low inclination of i = 95, where relativistic effects allow a complex interplay between the pencil and fan beam to produce the observed profiles. Although there is some systematic uncertainty in the derived inclination angle due to the complexity of the pulse profile modeling, we note that the inclination is in reasonable agreement with the face-on inclination of ≲8° inferred by studies of the orbit of the system that assume that the donor star is a 0.08 M_☉ hydrogen-depleted and partially degenerate star (Levine et al. 1988; Verbunt et al. 1990; Chakrabarty 1998).

The inclination is in moderate disagreement, however, with the $i\lesssim 33^\circ$ estimate for a 0.02 M_☉ helium or carbon–oxygen white dwarf donor (Verbunt et al. 1990; Chakrabarty 1998). This higher inclination case is supported by the presence of a complex of broad, double-peaked emission lines around 1 keV (Schulz et al. 2001; Krauss et al. 2007), which are consistent with an inclination in the range of 30°–40° (Schulz et al. 2013; Hemphill et al. 2017). The low inclination found by our pulse profile modeling can possibly be reconciled with the high inclination implied by the disk lines if the angular momenta of the accretion disk and neutron star are misaligned. This would result in a strong warp in the accretion disk, which could explain the disk flips that have been invoked to explain the torque reversals of 4U 1626−67 (van Kerkwijk et al. 1998; Wijers & Pringle 1999).

Due to the strong angular dependency of the cyclotron scattering cross section, we expect that the continuum emission and CRSF profile will depend on our viewing direction (see, e.g., Mészáros 1992). Here we perform a phase-resolved analysis to investigate how the spectral parameters change with phase. Phase-resolved spectroscopy of the NuSTAR data was explored in D17, however, their 20 phase bins had poor photon statistics, which precluded a detailed study of the fundamental CRSF. Motivated by our pulse profile modeling results, we instead divide the events into two phase intervals: ϕ = 0.0–0.5, which is more dominated by soft X-ray flux and contains the double peaks seen at low energies, and 0.5–1.0, where the profile is flat at low energies and more dominated by hard X-ray flux. Based on our pulse profile modeling, the phase 0.0–0.5 interval corresponds to viewing angles η of the first accretion column between 5° and 16°, while the 0.5–1.0 interval corresponds to angles between 17° and 22°.

3.4.1. Crab Ratio

Similar to our treatment of the phase-averaged data, we first performed a model-independent study of the phase-dependent data using Crab ratios. The normalized phase-resolved spectra by the canonical model of the Crab are shown in Figure 8. It is clear that the continuum emission and CRSF profile both depend strongly on the spin phase, with phases 0.5–1.0 being considerably harder spectrally than the 0.0–0.5 interval.

Zoom In Zoom Out Reset image size

3.4.2. Continuum and Cyclotron Line

We model the phase-resolved spectra with the modNPEX continuum models, which are successful in reproducing the phase-averaged spectra (Section 3.2.2). Due to the lower signal-to-noise ratio of the spectra, the iron line energy and width and those of the broad Gaussian component for the modNPEX model were fixed at the values from the phase-averaged spectrum (the NuSTAR results in Table 1).

Figure 9 shows the results of the spectral fitting of the spectrum accumulated for phase interval ϕ = 0.0–0.5. We do obtain a good fit with the modNPEX continuum and find that a single GABS component (Figure 9) is sufficient to reproduce the CRSF profile.

Zoom In Zoom Out Reset image size

We then fitted the spectrum from phase interval ϕ = 0.5–1.0 with the same model. The results are displayed in Figure 10. With a single GABS component, V-shaped residuals are still visible around the fundamental CRSF (see Figure 10(a)), similar to our results for the phase-averaged NuSTAR spectrum in Section 3.2.3. As before, we tried a CYAB model (a pseudo-Lorentzian profile) in place of the GABS component, but this produced worse fits compared with the GABS CRSF, with ${\chi }^{2}$/dof = 651.69/427 (Figure 10(b)). We then tried another GABS component (i.e., modNPEX × GABS × GABS) and an additive Gaussian component ([modNPEX + Gaussian] ×GABS). In both cases, the extra residuals around the CRSF were eliminated (Figures 10(c), (d)), and the Gaussian component normalization is negative at the 5σ level. The probability of this feature arising by chance, determined by simulating 400,000 data sets with simftest in XSPEC, is 1.8 × 10⁻⁵. The best-fit parameters using the modNPEX model are summarized in Table 4.

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Parameters for the Phase 0.0–0.5 and 0.5–1.0 NuSTAR Spectra of 4U 1626−67

	modNPEX	modNPEX	modNPEX	(modNPEX + Gaussian)
	× GABS	× GABS	× GABS × GABS	× GABS
Pulse phase ϕ	0.0–0.5	0.5–1.0
kTBB (keV)	${0.49}_{-0.04}^{+0.04}$	${0.42}_{-0.07}^{+0.07}$	${0.42}_{-0.07}^{+0.07}$	${0.42}_{-0.07}^{+0.07}$
ABBa	${244}_{-85}^{+158}$	${352}_{-214}^{+1165}$	${359}_{-210}^{+1026}$	${367}_{-217}^{+1079}$
α	${0.47}_{-0.01}^{+0.01}$	${1.02}_{-0.09}^{+0.08}$	${1.00}_{-0.10}^{+0.08}$	${1.02}_{-0.09}^{+0.08}$
kTNPEX (keV)	${5.6}_{-0.2}^{+0.3}$	${6.0}_{-0.1}^{+0.1}$	${6.0}_{-0.1}^{+0.1}$	${6.0}_{-0.1}^{+0.1}$
An(×10−2)b	${1.9}_{-0.2}^{+0.2}$	${4.8}_{-0.7}^{+0.6}$	${4.7}_{-0.7}^{+0.7}$	${4.9}_{-0.7}^{+0.7}$
Ap(×10−5)b	${3.2}_{-1.2}^{+1.1}$	${10.3}_{-0.3}^{+0.3}$	${12.6}_{-0.3}^{+0.3}$	${10.4}_{-0.3}^{+0.3}$
${E}_{\mathrm{Gaussian}}$ (keV)	19.8(fix)	19.8(fix)	19.8(fix)	19.8(fix)
σGaussian (keV)	4.8(fix)	4.8(fix)	4.8(fix)	4.8(fix)
AGaussian(×10−3)	${1.1}_{-0.3}^{+0.3}$	${5.1}_{-0.3}^{+0.3}$	${5.5}_{-0.3}^{+0.3}$	${5.4}_{-0.3}^{+0.3}$
EFe (keV)	6.76(fix)	6.76(fix)	6.76(fix)	6.76(fix)
σFe (keV)	0.15(fix)	0.15(fix)	0.15(fix)	0.15(fix)
AFe(×10−4)c	${1.2}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$
Ea (keV)	${36.5}_{-0.4}^{+0.5}$	${37.9}_{-0.1}^{+0.1}$	${39.2}_{-0.7}^{+1.0}$	${38.9}_{-0.5}^{+0.5}$
σCRSF (keV)	${4.6}_{-0.4}^{+0.4}$	${4.0}_{-0.1}^{+0.1}$	${2.8}_{-0.5}^{+0.5}$	${3.0}_{-0.3}^{+0.3}$
dCRSF	${10.8}_{-1.6}^{+1.7}$	${15.9}_{-0.7}^{+0.7}$	${10.3}_{-5.5}^{+2.8}$	${11.7}_{-2.1}^{+1.7}$
Eabs (keV)	⋯	⋯	${34.2}_{-1.5}^{+2.2}$	${32.5}_{-0.8}^{+0.7}$
σabs (keV)	⋯	⋯	${2.9}_{-0.5}^{+0.8}$	${2.8}_{-0.5}^{+0.5}$
dabs	⋯	⋯	${4.8}_{-2.6}^{+5.6}$	⋯
Normabs × 10−3	⋯	⋯	⋯	$-{1.4}_{-0.5}^{+0.5}$
χ2/dof	453.52/427	522.50/427	493.80/424	494.55.47/424
pnulld	1.8 × 10−1	1.1 × 10−3	1.1 × 10−2	1.0 × 10−2

Notes.

^a ${R}_{\mathrm{km}}^{2}/{d}_{10}^{2}$, where R_km is the radius of the blackbody in kilometers and d₁₀ is the distance in units of 10 kpc. ^bNormalization of negative (n) and positive (p) NPEX components, defined at 1 keV in units of photons keV⁻¹ cm⁻² s⁻¹. ^cNormalization of the Gaussian, defined in units of photons keV⁻¹ cm⁻² s⁻¹. ^dNull hypothesis probability.

Download table as:  ASCIITypeset image

Hence, we conclude that we have again tentatively detected the distorted CRSF, which cannot be represented by a simple Gaussian or pseudo-Lorentzian absorption model, similar to the CRSF seen in the phase-averaged NuSTAR spectrum.

3.4.3. Dim Phase Spectrum

Iwakiri et al. (2012) reported the possible detection of an emission line-like feature in the dim phase (ϕ = 0.125–0.250 in Figure 5) of 4U 1626−67 during the spin-down state. We extracted a spectrum from this phase interval and fitted it with the modNPEX model; the results are displayed in Figure 11(a) and Table 5. The spectrum is well reproduced by only the continuum model—the CRSF feature does not appear in the dim phase significantly. However, the residuals around 45 keV do indicate a possible emission line-like feature, although the signal is weak compared with the background level. If we add a Gaussian emission component in the same manner as Iwakiri et al. (2012), with the center energy fixed at their best-fit value E = 40.5 keV (Figure 11(b)), the fitted width of $\sigma \,={8.5}_{-2.5}^{+3.8}\,\mathrm{keV}$ is consistent with the width found in the spin-down state, which showed σ > 4.5 keV, whereas the normalization of ${2.2}_{-0.7}^{+0.8}\times {10}^{-4}$ photons cm⁻² s⁻¹ is lower than the spin-down observation's ${8.8}_{-2.9}^{+7.2}\times {10}^{-4}$ ph cm⁻² s⁻¹. The chance probability of improvements of adding this emission feature is 2.0 × 10⁻², determined by simulating 20,000 data sets with simftest in XSPEC. Therefore, we concluded that we have only marginally detected the emission line at about the 2σ level. The observed unabsorbed fluxes for the dim phase of the spin-down state and the spin-up state in the 3–60 keV band are ${2.2}_{-0.2}^{+0.2}\times {10}^{-10}$ and ${6.0}_{-1.0}^{+0.1}\times {10}^{-10}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, respectively (when we calculate the fluxes, the blackbody parameters are fixed). We note that while D17 claimed the presence of an absorption feature in the dim phase of the pulse, their chosen phase interval is different from ours, corresponding to approximately ϕ = 0.2–0.35 in Figure 5, and thus this discrepancy is not problematic.

Zoom In Zoom Out Reset image size

Table 5.  Best-fit Parameters for the Phase 0.125–0.250 NuSTAR Spectra of 4U 1626–67

	modNPEX	modNPEX
		+ Gaussian emission
kTBB (keV)	${0.41}_{-0.05}^{+0.05}$	${0.47}_{-0.05}^{+0.05}$
ABBa	${684}_{-397}^{+1250}$	${324}_{-159}^{+424}$
α	${0.10}_{-0.09}^{+0.10}$	$-{0.15}_{-0.04}^{+0.18}$
kTNPEX (keV)	${5.3}_{-0.5}^{+0.5}$	${4.6}_{-0.6}^{+6.0}$
An(×10−2)b	${1.6}_{-0.2}^{+0.2}$	${1.3}_{-0.2}^{+0.3}$
Ap(×10−5)b	${0.5}_{-0.4}^{+1.1}$	<4.9
EGaussian (keV)	19.8(fix)	19.8(fix)
σGaussian (keV)	4.8(fix)	4.8(fix)
AGaussian(×10−3)	${0.5}_{-0.2}^{+0.2}$	${1.0}_{-0.3}^{+0.3}$
EFe (keV)	6.76(fix)	6.76(fix)
σFe (keV)	0.15(fix)	0.15(fix)
AFe(×10−4)c	${0.8}_{-0.5}^{+0.5}$	${0.9}_{-0.5}^{+0.5}$
Ea (keV)	⋯	40.5(fix)
σ (keV)	⋯	${8.5}_{-2.5}^{+3.8}$
Norm (×10−4)	⋯	${2.2}_{-0.7}^{+0.7}$
χ2/dof	395.68/384	366.99/382
pnulld	3.3 × 10−1	7.0 × 10−1

Notes.

^a ${R}_{\mathrm{km}}^{2}/{d}_{10}^{2}$, where R_km is the radius of the blackbody in kilometers and d₁₀ is the distance in units of 10 kpc. ^bNormalization of negative (n) and positive (p) NPEX components, defined at 1 keV in units of photons keV⁻¹ cm⁻² s⁻¹. ^cNormalization of the Gaussian, defined in units of photons keV⁻¹ cm⁻² s⁻¹. ^dNull hypothesis probability.

Download table as:  ASCIITypeset image

We have found that the continuum emission is different between the 2010 Suzaku and 2015 NuSTAR observations. Let us consider the implications of this result.

According to Ghosh & Lamb (1979), for accretion via a disk, the rate of change in the pulse period, $\dot{P}$, is proportional to ${\dot{M}}^{6/7}$, where $\dot{M}$ is the total mass accretion rate onto the neutron star. Thus, we can estimate the accretion rate onto 4U 1626−67 using $\dot{P}$ from Fermi/GBM monitoring. As found in Section 2.3, $\dot{P}$ over the Suzaku observation was −2.8 × 10⁻¹¹ s s⁻¹, while over the NuSTAR observation it was −3.3 × 10⁻¹¹ s s⁻¹. The 15% higher $\dot{P}$ during the NuSTAR observation would thus imply a 13% higher accretion rate per Ghosh & Lamb (1979).

Based on the Crab ratios of these data sets (Figure 2), as well as our spectral fitting results (Table 1), hard X-ray photons are suppressed in the NuSTAR observation relative to the 2010 Suzaku spectrum. Because hard X-ray photons are mainly produced by thermal comptonization in the accretion column (Becker & Wolff 2007), we can infer that the electron temperature of the plasma decreased. However, the photon index α did not change significantly between the observations.

Summarizing our interpretations of our timing and spectral results, the NuSTAR observation saw a higher accretion rate and lower temperature in the accretion column compared with the Suzaku observation. At least in terms of these results, the spectral variation between the 2010 Suzaku and 2015 NuSTAR observations indicates that the decrease in the plasma temperature is caused by the increased accretion rate. However, the relation between the accretion rate and the temperature is not simple due to the complex of radiative transfer processes under the strong magnetic field. Future theoretical study is needed to verify our interpretation.

The high sensitivity of NuSTAR in the hard X-ray band and our analyses of the phase-averaged and phase-resolved spectra lead us to conclude that the most appropriate empirical model for the continuum emission of 4U 1626−67 is the modNPEX model. The results indicate that the extra Gaussian is needed to allow the NPEX model-based fit to provide a good description of the physical continuum. Similar broad Gaussian features have been detected from several other accretion-powered pulsars, e.g., 4U 0115+63, Cen X-3, Her X-1, and A 0535+26 (Suchy et al. 2008; Ferrigno et al. 2009; Vasco et al. 2013; Ballhausen et al. 2017). Because the broad Gaussian feature appears around the spectral peak, it is inferred that the main difference between the simple NPEX continuum model and the observed spectrum is the spectral shape of the quasi-exponential cutoff. The NPEX model approximates the Wien peak-like cutoff, which is expected for a pure thermal comptonization regime, but the shape of the spectral cutoff is mainly determined by the contribution ratio of the bulk and thermal comptonization in their photon propagation process. Moreover, because there are three types of seed photons (bremsstrahlung and cyclotron emission created along the column, blackbody emission from the base of the column), the cutoff shape is also related to the contribution ratio of the seed photons. Therefore, a physical model fitting is important to investigate the origin of the broad Gaussian feature observed from 4U 1626−67. However, this is beyond the scope of the present paper. While D17 did fit 4U 1626–67 with the Becker & Wolff (2007) bulk and thermal comptonization model, their description includes an additional reflection component and the unusual second harmonic absorption feature described earlier. Because of this and because the individual spectral contributions of the different types of comptonized seed photons are not shown, we cannot interpret the 20 keV residual within the D17 picture.

Using our new relativistic ray tracing code, we successfully reproduce the energy-resolved pulse profiles obtained by NuSTAR (Section 3.3). We find that the beam pattern is energy dependent. The qualitative explanation for the observed pulse profile evolution is found in the anisotropy of the Thomson scattering cross section in a strongly magnetized plasma. In this case the cross section depends on the photon energy and is different for the ordinary and extraordinary polarization modes. Kii et al. (1986) simulated the 4U 1626−67 pulse profile during spin-up, as obtained by Tenma during 1983 May 3–5. Their simulation is based on anisotropic radiation transfer calculations by Nagel (1981), ignoring vacuum polarization. For a cylindrical geometry with a diameter of the Thomson scattering optical depth τ_Th = 80, Kii et al. (1986) find the following: first, ordinary-mode photons dominate in the lower energy band, whereas extraordinary-mode photons become dominant toward higher energies. Second, the maximum emission angle for extraordinary-mode photons is offset with respect to the magnetic field, and that offset becomes slightly stronger at higher energies. Third, the emission profile with respect to the magnetic field becomes wider at higher energies. Comparing these with our pulse profile modeling results, we can thus draw qualitative connections between our fan-beam evolution and Kii et al.'s (1986) ordinary-mode emission profile and between our pencil-beam evolution and their extraordinary-mode emission profile. Therefore, we suggest that the origin of the observed pulse profile evolution is found in the energy and polarization dependence of the scattering cross section in a strongly magnetized plasma.

X-ray polarimetry will be useful to validate this hypothesis, as ordinary- and extraordinary-mode photons contribute to polarization with opposite signs. Consider, as an example, the future mission IXPE (Weisskopf et al. 2016), which will operate in the 2–8 keV band. Under our suggestion that the fan and pencil beams, respectively, correspond to ordinary- and extraordinary-mode photons, Figure 5 predicts that in IXPE's band, the polarization degree will be nearly zero around phase ϕ = 0.45, where the pencil and fan contributions are comparable, and maximal around phases ϕ = 0.3 and ϕ = 0.75, where one beam dominates over the other. Meanwhile, the hard X-ray band covered by, e.g., X-Calibur (Beilicke et al. 2014) should find higher overall polarization with a minimum around phase ϕ = 0.7.

We have found that the observed fundamental CRSF in the phase-averaged and the phase ϕ = 0.5–1.0 spectra is better described with a two-Gaussian absorption structure than with a single Gaussian or pseudo-Lorentzian profile (Section 3.4.2). To highlight this asymmetry, we show in Figure 12 the ratio between the data and the continuum component of the best-fit model, where we have excluded the CRSF model. The best-fit two-Gaussian absorption model is overlaid, showing broadening at lower energies. This is similar to what was reported in the NuSTAR observation of Cep X-4 (Fürst et al. 2015). Some explanations of this distorted profile have been proposed. For example, Nishimura (2011) proposed that the fundamental line profile becomes asymmetric and shallower toward lower energies due to the superposition of multiple CRSFs produced at different altitudes along the column, assuming some gradient in the density, temperature, and magnetic field. Another interpretation is photon spawning due to inelastic scattering at higher harmonics (Schönherr et al. 2007)—as electrons excited into higher Landau levels cascade down to the ground state, many emitted photons will have similar energy to the fundamental, "filling in" the fundamental line. In addition, our pulse profile modeling suggests that in the phase interval ϕ = 0.5–1.0, the second accretion column contributes about 10% of the flux of the first column. Thus, the emission from the second pole may slightly contribute to the CRSF shape. However, a detailed study of the spectra resulting from the mixing of light from the two columns is beyond the scope of this work.

Zoom In Zoom Out Reset image size

Our phase-resolved spectral analysis also suggests that the CRSF profile depends on the spin phase. We have found that the CRSF in the phase ϕ = 0.0–0.5 spectra, which corresponds to the angle η_AC1 = 5°–16° (Figure 7), is well reproduced by a simple GABS model. Figure 13 shows a comparison of the CRSF shape between phase ϕ = 0.5–1.0 and ϕ = 0.0–0.5 in terms of data-to-model ratio (with the absorption components excluded). The depth of the line at ϕ = 0.0–0.5 is clearly shallower than at ϕ = 0.5–1.0. It is noteworthy that the wings of the two profiles resemble each other, while the core of the ϕ = 0.5–1.0 profile appears to be filled in compared with the ϕ = 0.0–0.5 profile. To evaluate the width of the CRSF profiles, we calculate the FWHM from these ratios. The results are 12.9 and 11.4 keV, corresponding to ϕ = 0.0–0.5 and ϕ = 0.5–1.0 , respectively. Therefore, the CRSF is wider and shallower between η_AC1 = 5°–16° compared with η_AC1 = 17°–21° (see Figure 7). The relation between the CRSF shape and the angle η is qualitatively consistent with theoretical simulations when assuming a slab-type geometry (Isenberg et al. 1998; Schwarm et al. 2017a). The relation is also qualitatively consistent with the theoretical results assuming a cylindrical geometry illuminated by anisotropic injections and magnetic field gradients (Nishimura 2015).

Zoom In Zoom Out Reset image size

The absorption feature disappears in the dim phase spectrum (phase ϕ = 0.125–0.250). Although an emission line-like feature was detected in the Suzaku spectrum observed during the 2006 spin-down state of 4U 1626−67 (Iwakiri et al. 2012), it is only marginally detected in the NuSTAR observation (the probability of the feature arising by chance is 2.0 × 10⁻²). This fact indicates that the flux of the emission component decreased to the point of nondetectability in the 2015 spin-up state. According to our spectral fits (see Section 3.4.3), the intensity of the emission feature during the NuSTAR observation is about four times lower than during the 2006 Suzaku observation. On the other hand, the 3–60 keV flux of the NuSTAR observation is about three times higher than that of the 2006 Suzaku observation. As a consistency check we refitted the 2006 Suzaku observation with the same modNPEX model used for the NuSTAR observation in Section 3.4.3 (Figure 14). We find consistent results with the emission line-like feature around 40 keV being possibly detected (the probability of this feature arising by chance is 8.0 × 10⁻⁴ using 20,000 data sets with simftest) and the fitted continuum approximating the exponentially cut-off power law used by Iwakiri et al. (2012) because the contribution of the broad Gaussian and the positive exponential component are comparatively small (A_Gaussian < 0.4 × 10⁻³ keV⁻¹ cm⁻² s⁻¹ and A_p <1.5 × 10⁻⁵ keV⁻¹ cm⁻² s⁻¹). According to Section 3.3 the dim phase corresponds to η_AC1 ≈ 5°, a viewing angle almost parallel to the magnetic field. Nishimura (2015) specifically addressed 4U 1626−67 and found in simulations that an emission feature could arise at around 50 keV for viewing angles nearly parallel to the magnetic field, which is qualitatively consistent with our observations as well.

Zoom In Zoom Out Reset image size