Spectral and Timing Analysis of the accretion-powered pulsar 4U 1626-67 observed with Suzaku and NuSTAR

TL;DR

This study analyzes spectral and timing data of the pulsar 4U 1626-67 from Suzaku and NuSTAR, revealing spectral changes linked to accretion rate variations and providing insights into the geometry of accretion columns through relativistic ray tracing.

Contribution

It introduces a combined spectral and pulse profile analysis using relativistic ray tracing to infer accretion column geometry in a pulsar, linking spectral changes to accretion rate variations.

Findings

Spectral shape changed between observations, with flux increase in 3-10 keV and decrease in 30-60 keV.

Detected asymmetry in the fundamental cyclotron line profile.

Derived geometrical parameters of accretion columns from pulse profile decomposition.

Abstract

We present an analysis of the spectral shape and pulse profile of the accretion-powered pulsar 4U 1626-67 observed with Suzaku and NuSTAR during a spin-up state. The pulsar, which experienced a torque reversal to spin-up in 2008, has a spin period of 7.7 s. Comparing the phase-averaged spectra obtained with Suzaku in 2010 and with NuSTAR in 2015, we find that the spectral shape changed between the two observations: the 3-10 keV flux increased by 5% while the 30-60 keV flux decreased significantly by 35%. Phase-averaged and phase-resolved spectral analysis shows that the continuum spectrum observed by NuSTAR is well described by an empirical NPEX continuum with an added broad Gaussian emission component around the spectral peak at 20 keV. Taken together with the observed Pdot value obtained from Fermi/GBM, we conclude that the spectral change between the Suzaku and NuSTAR observations was likely caused by an increase of the accretion rate. We also report the possible detection of asymmetry in the profile of the fundamental cyclotron line. Furthermore, we present a study of the energy-resolved pulse profiles using a new relativistic ray tracing code, where we perform a simultaneous fit to the pulse profiles assuming a two-column geometry with a mixed pencil- and fan-beam emission pattern. The resulting pulse profile decompositions enable us to obtain geometrical parameters of accretion columns (inclination, azimuthal and polar angles) and a fiducial set of beam patterns. This information is important to validate the theoretical predictions from radiation transfer in a strong magnetic field.

Tables

Table 1. 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 $\times$ GABS
	Suzaku	NuSTAR
$k{T}_{\mathrm{BB}}$ (keV)	${0.47}_{-0.06}^{+0.06}$	${0.47}_{-0.04}^{+0.04}$
$A_{\mathrm{BB}}$aa$R_{\mathrm{km}}^2/d_{10}^2$, where $R_{\mathrm{km}}$ is the radius of the black-body in km and $d_{10}$ is the distance in units of 10 kpc.	${258}_{-104}^{+271}$	${254}_{-85}^{+174}$
$\alpha$	${0.53}_{-0.18}^{+0.16}$	${0.47}_{-0.10}^{+0.09}$
$k{T}_{\mathrm{NPEX}}$ (keV)	${7.1}_{-0.4}^{+0.6}$	${5.9}_{-0.1}^{+0.1}$
$A_n(\times {10}^{-2})$bbNormalization of the power-law. Defined at 1 keV in units of photons keV-1 cm-2 s-1.	${2.3}_{-0.6}^{+0.7}$	${2.6}_{-0.4}^{+0.4}$
$A_p(\times {10}^{-5})$bbNormalization of the power-law. Defined at 1 keV in units of photons keV-1 cm-2 s-1.	${4.0}_{-1.8}^{+1.4}$	${7.0}_{-1.3}^{+0.9}$
$E_{\mathrm{Gaussian}}$ (keV)	19.8(fix)	${19.8}_{-0.5}^{+0.4}$
${\sigma }_{\mathrm{Gaussian}}$ (keV)	4.8(fix)	${4.8}_{-0.5}^{+0.6}$
$A_{\mathrm{Gaussian}}(\times {10}^{-3})$ddA cross calibration constant for XIS-BI relative to XIS-FI.	${1.8}_{-1.7}^{+2.3}$	${3.1}_{-0.8}^{+1.2}$
$E_{\mathrm{Fe}}$ (keV)	${6.75}_{-0.06}^{+0.06}$	${6.76}_{-0.05}^{+0.05}$
${\sigma }_{\mathrm{Fe}}$ (keV)	${0.13}_{-0.08}^{+0.06}$	${0.15}_{-0.06}^{+0.06}$
$A_{\mathrm{Fe}}(\times {10}^{-4})$ccNormalization of the Gaussian. Defined in units of photons keV-1 cm-2 s-1.	${1.3}_{-0.4}^{+0.5}$	${1.2}_{-0.2}^{+0.3}$
$E_{\mathrm{CRSF}}$ (keV)	${38.2}_{-0.8}^{+0.9}$	${37.7}_{-0.1}^{+0.1}$
${\sigma }_{\mathrm{CRSF}}$ (keV)	${5.1}_{-0.7}^{+0.7}$	${4.2}_{-0.1}^{+0.1}$
$d_{\mathrm{CRSF}}$	${22.9}_{-2.6}^{+5.7}$	${14.7}_{-0.8}^{+0.8}$
$C_{{\mathrm{XIS}}_{\mathrm{BI}}}$ddA cross calibration constant for XIS-BI relative to XIS-FI.	0.941${}_{-0.007}{}^{}{}_{}{}^{+0.007}$	—
$C_{\mathrm{FPMB}}$eeA cross calibration constant for FPMB relative to FPMA.	—	1.002${}_{-0.002}{}^{}{}_{}{}^{+0.002}$
${\chi }^2$ / d.o.f	323.48/316	533.08/422
$p_{\mathrm{null}}$ffNull hypothesis probability.	3.7$\times$10-1	1.9$\times$10-4

Table 2. 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)
	$\times$ CYAB	$\times$ GABS $\times$ GABS	$\times$ GABS
$k{T}_{\mathrm{BB}}$ (keV)	${0.47}_{-0.04}^{+0.04}$	${0.47}_{-0.03}^{+0.03}$	${0.47}_{-0.03}^{+0.03}$
$A_{\mathrm{BB}}$aa$R_{\mathrm{km}}^2/d_{10}^2$, where $R_{\mathrm{km}}$ is the radius of the black-body in km and $d_{10}$ is the distance in units of 10 kpc.	${263}_{-89}^{+167}$	${266}_{-86}^{+154}$	${266}_{-87}^{+155}$
$\alpha$	${0.52}_{-0.07}^{+0.06}$	${0.48}_{-0.05}^{+0.05}$	${0.48}_{-0.05}^{+0.05}$
$k{T}_{\mathrm{NPEX}}$ (keV)	${6.3}_{-0.1}^{+0.1}$	${5.9}_{-0.1}^{+0.1}$	${5.9}_{-0.1}^{+0.1}$
$A_n(\times {10}^{-2})$bbNormalization of the power-law. Defined at 1 keV in units of photons keV-1 cm-2 s-1.	${2.6}_{-0.4}^{+0.3}$	${2.6}_{-0.2}^{+0.2}$	${2.6}_{-0.2}^{+0.2}$
$A_p(\times {10}^{-5})$bbNormalization of the power-law. Defined at 1 keV in units of photons keV-1 cm-2 s-1.	${6.1}_{-1.0}^{+0.7}$	${7.1}_{-0.4}^{+0.4}$	${7.1}_{-0.4}^{+0.4}$
$E_{\mathrm{Gaussian}}$ (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)
$A_{\mathrm{Gaussian}}(\times {10}^{-3})$ccNormalization of the Gaussian. Defined in units of photons keV-1 cm-2 s-1.	${3.1}_{-0.8}^{+1.3}$	3.1${}_{-0.1}{}^{}{}_{}{}^{+0.3}$	3.1${}_{-0.1}{}^{}{}_{}{}^{+0.3}$
$E_{\mathrm{Fe}}$ (keV)	${6.76}_{-0.05}^{+0.05}$	${6.76}_{-0.05}^{+0.05}$	${6.76}_{-0.05}^{+0.05}$
${\sigma }_{\mathrm{Fe}}$ (keV)	${0.15}_{-0.06}^{+0.06}$	${0.15}_{-0.06}^{+0.06}$	${0.15}_{-0.06}^{+0.06}$
$A_{\mathrm{Fe}}(\times {10}^{-4})$ccNormalization of the Gaussian. Defined in units of photons keV-1 cm-2 s-1.	${1.2}_{-0.2}^{+0.3}$	${1.2}_{-0.2}^{+0.3}$	${1.2}_{-0.2}^{+0.3}$
$E_{\mathrm{CRSF}}$ (keV)	${36.6}_{-0.1}^{+0.1}$	${39.9}_{-1.3}^{+0.6}$	${39.2}_{-0.7}^{+0.8}$
${\sigma }_{\mathrm{CRSF}}$ (keV)	${5.3}_{-0.3}^{+0.3}$	${2.6}_{-0.8}^{+0.6}$	${2.9}_{-0.5}^{+0.5}$
${\tau }_{\mathrm{CRSF}}$	${1.6}_{-0.1}^{+0.1}$	—	—
$d_{\mathrm{CRSF}}$	—	${5.3}_{-3.8}^{+5.9}$	${8.8}_{-3.2}^{+2.8}$
$E_{\mathrm{abs}}$ (keV)	—	${35.9}_{-2.7}^{+1.4}$	${32.8}_{-1.1}^{+0.9}$
${\sigma }_{\mathrm{abs}}$ (keV)	—	${3.6}_{-0.7}^{+0.5}$	${3.4}_{-0.5}^{+0.5}$
$d_{\mathrm{abs}}$	—	${8.7}_{-5.6}^{+4.2}$	—
$A_{\mathrm{abs}}(\times {10}^{-3})$ccNormalization of the Gaussian. Defined in units of photons keV-1 cm-2 s-1.	—	—	$-{1.2}_{-0.4}^{+0.4}$
$C_{\mathrm{FPMB}}$ddA cross calibration constant for FPMB relative to FPMA.	1.002${}_{-0.002}{}^{}{}_{}{}^{+0.002}$	1.002${}_{-0.002}{}^{}{}_{}{}^{+0.002}$	1.002${}_{-0.002}{}^{}{}_{}{}^{+0.002}$
${\chi }^2$ / d.o.f	596.91/422	508.08/421	507.98/421
$p_{\mathrm{null}}$eeNull hypothesis probability.	2.1$\times$10-9	2.3$\times$10-3	2.3$\times$10-3

Table 3. Table 3: Fit parameters of the accretion column model. Parameters not listed here are fixed or tied. In particular we impose r AC1 = r AC2 subscript 𝑟 AC1 subscript 𝑟 AC2 r_{\mathrm{AC1}}=r_{\mathrm{AC2}} and h AC1 = h AC2 subscript ℎ AC1 subscript ℎ AC2 h_{\mathrm{AC1}}=h_{\mathrm{AC2}} for the column size, N p 2 = N p 1 subscript 𝑁 subscript p 2 subscript 𝑁 subscript p 1 N_{\mathrm{p}_{2}}=N_{\mathrm{p}_{1}} , σ p 2 = σ p 1 subscript 𝜎 subscript p 2 subscript 𝜎 subscript p 1 \sigma_{\mathrm{p}_{2}}=\sigma_{\mathrm{p}_{1}} , η ¯ p 2 = η ¯ p 1 subscript ¯ 𝜂 subscript p 2 subscript ¯ 𝜂 subscript p 1 \bar{\eta}_{\mathrm{p}_{2}}=\bar{\eta}_{\mathrm{p}_{1}} for the pencil beam, and σ f 2 = σ f 1 subscript 𝜎 subscript f 2 subscript 𝜎 subscript f 1 \sigma_{\mathrm{f}_{2}}=\sigma_{\mathrm{f}_{1}} and η ¯ f 2 = η ¯ f 1 = 90 ∘ subscript ¯ 𝜂 subscript f 2 subscript ¯ 𝜂 subscript f 1 superscript 90 \bar{\eta}_{\mathrm{f}_{2}}=\bar{\eta}_{\mathrm{f}_{1}}=90^{\circ} for the fan beam. The given errors correspond to the 90% confidence level.

global	column1					column2
$i$	${\Phi }_{\mathrm{AC1}}$	${\Theta }_{\mathrm{AC1}}$	$r_{\mathrm{AC1}}$	$h_{\mathrm{AC1}}$		${\Phi }_{\mathrm{AC2}}$	${\Theta }_{\mathrm{AC2}}$	${\chi }^2$/d.o.f.
[deg]	[deg]	[deg]	[m]	[m]		[deg]	[deg]
$9.5\pm 0.2$	$73.6\pm 0.3$	${12.50}_{-0.11}^{+0.07}$	$590\pm 50$	${252}_{-14}^{+12}$		$290\pm 2$	$158.1\pm 0.2$	$1146/631$
	pencil1			fan1		fan2
$\Delta E$	$N_{p_1}$	${\sigma }_{p_1}$	${\overline{\eta }}_{p_1}$	$N_{f_1}$	${\sigma }_{f_1}$	$N_{f_2}$
[keV]	[cts/s]	[deg]	[deg]	[cts/s]	[deg]	[cts/s]
$3$ – $5$	${0.95}_{-0.02}^{+0.01}$	${4.21}_{-0.06}^{+0.07}$	$5.73\pm 0.06$	${1.21}_{-0.02}^{+0.03}$	${180}_{-16}^{+0}$	${0.011}_{-0.011}^{+0.013}$
$5$ – $7$	$1.14\pm 0.02$	$4.37\pm 0.06$	$6.12\pm 0.06$	${1.09}_{-0.02}^{+0.03}$	${180}_{-26}^{+0}$	$0.08\pm 0.02$
$7$ – $9$	${1.03}_{-0.02}^{+0.03}$	$4.66\pm 0.05$	${6.65}_{-0.06}^{+0.07}$	${0.84}_{-0.04}^{+0.03}$	${180}_{-67}^{+0}$	$0.16\pm 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\pm 0.02$
$11$ – $13$	${0.77}_{-0.09}^{+0.04}$	$7.3\pm 0.7$	$9.4\pm 0.2$	${0.17}_{-0.04}^{+0.08}$	${20}_{-6}^{+12}$	$0.13\pm 0.02$
$13$ – $15$	${0.56}_{-0.03}^{+0.02}$	$7.8\pm 0.6$	${10.7}_{-0.2}^{+0.3}$	${0.11}_{-0.02}^{+0.03}$	${14.8}_{-0.7}^{+4.0}$	$0.11\pm 0.02$
$15$ – $17$	${0.46}_{-0.02}^{+0.03}$	$8.4\pm 0.5$	$12.4\pm 0.4$	$0.06\pm 0.02$	${13.0}_{-1.7}^{+0.9}$	$0.082\pm 0.009$
$17$ – $20$	${0.55}_{-0.04}^{+0.02}$	$8.3\pm 0.3$	${14.8}_{-0.6}^{+0.3}$	${0.026}_{-0.008}^{+0.050}$	$7_{-3}^{+6}$	$0.08\pm 0.01$
$20$ – $25$	$0.53\pm 0.01$	$7.8\pm 0.2$	$16.2\pm 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\pm 0.3$	${17.0}_{-0.7}^{+0.5}$	${0.02}_{-0.01}^{+0.02}$	${11}_{-5}^{+2}$	$0.025\pm 0.005$
$30$ – $50$	${0.070}_{-0.008}^{+0.004}$	${8.3}_{-0.5}^{+0.6}$	$18.1\pm 0.7$	$\le 0.007$	${14}_{-8}^{+167}$	$0.018\pm 0.004$

Table 4. 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)
	$\times$ GABS	$\times$ GABS	$\times$ GABS $\times$ GABS	$\times$ GABS
Pulse phase $\varphi$	0.0–0.5	0.5–1.0
$k{T}_{\mathrm{BB}}$ (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}$
$A_{\mathrm{BB}}$aa$R_{\mathrm{km}}^2/d_{10}^2$, where $R_{\mathrm{km}}$ is the radius of the black-body in km and $d_{10}$ is the distance in units of 10 kpc.	${244}_{-85}^{+158}$	${352}_{-214}^{+1165}$	${359}_{-210}^{+1026}$	${367}_{-217}^{+1079}$
$\alpha$	${0.47}_{-0.01}^{+0.01}$	${1.02}_{-0.09}^{+0.08}$	${1.00}_{-0.10}^{+0.08}$	${1.02}_{-0.09}^{+0.08}$
$k{T}_{\mathrm{NPEX}}$ (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}$
$A_n(\times {10}^{-2})$bbNormalization of negative ($n$) and positive ($p$) NPEX components, defined at 1 keV in units of photons keV-1 cm-2 s-1.	${1.9}_{-0.2}^{+0.2}$	${4.8}_{-0.7}^{+0.6}$	${4.7}_{-0.7}^{+0.7}$	${4.9}_{-0.7}^{+0.7}$
$A_p(\times {10}^{-5})$bbNormalization of negative ($n$) and positive ($p$) NPEX components, defined at 1 keV in units of photons keV-1 cm-2 s-1.	${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)
${\sigma }_{\mathrm{Gaussian}}$ (keV)	4.8(fix)	4.8(fix)	4.8(fix)	4.8(fix)
$A_{\mathrm{Gaussian}}(\times {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}$
$E_{Fe}$ (keV)	6.76(fix)	6.76(fix)	6.76(fix)	6.76(fix)
${\sigma }_{Fe}$ (keV)	0.15(fix)	0.15(fix)	0.15(fix)	0.15(fix)
$A_{\mathrm{Fe}}(\times {10}^{-4})$ccNormalization of the Gaussian. Defined in units of photons keV-1 cm-2 s-1.	${1.2}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$
$E_a$ (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}$
${\sigma }_{\mathrm{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}$
$d_{\mathrm{CRSF}}$	${10.8}_{-1.6}^{+1.7}$	${15.9}_{-0.7}^{+0.7}$	${10.3}_{-5.5}^{+2.8}$	${11.7}_{-2.1}^{+1.7}$
$E_{\mathrm{abs}}$ (keV)	—	—	${34.2}_{-1.5}^{+2.2}$	${32.5}_{-0.8}^{+0.7}$
${\sigma }_{\mathrm{abs}}$ (keV)	—	—	${2.9}_{-0.5}^{+0.8}$	${2.8}_{-0.5}^{+0.5}$
$d_{\mathrm{abs}}$	—	—	${4.8}_{-2.6}^{+5.6}$	—
Norm${}_{abs}{}^{}\times {10}^{-3}$	—	—	—	$-{1.4}_{-0.5}^{+0.5}$
${\chi }^2$ / d.o.f	453.52/427	522.50/427	493.80/424	494.55.47/424
$p_{\mathrm{null}}$ddNull hypothesis probability.	1.8$\times$10-1	1.1$\times$10-3	1.1$\times$10-2	1.0$\times$10-2

Table 5. Table 5: Best fit parameters for the phase 0.125–0.250 NuSTAR spectra of 4U 1626-67.

	modNPEX	modNPEX
		+ Gaussian emission
$k{T}_{\mathrm{BB}}$ (keV)	${0.41}_{-0.05}^{+0.05}$	${0.47}_{-0.05}^{+0.05}$
$A_{\mathrm{BB}}$aa$R_{\mathrm{km}}^2/d_{10}^2$, where $R_{\mathrm{km}}$ is the radius of the black-body in km and $d_{10}$ is the distance in units of 10 kpc.	${684}_{-397}^{+1250}$	${324}_{-159}^{+424}$
$\alpha$	${0.10}_{-0.09}^{+0.10}$	$-{0.15}_{-0.04}^{+0.18}$
$k{T}_{\mathrm{NPEX}}$ (keV)	${5.3}_{-0.5}^{+0.5}$	${4.6}_{-0.6}^{+6.0}$
$A_n(\times {10}^{-2})$bbNormalization of negative ($n$) and positive ($p$) NPEX components, defined at 1 keV in units of photons keV-1 cm-2 s-1.	${1.6}_{-0.2}^{+0.2}$	${1.3}_{-0.2}^{+0.3}$
$A_p(\times {10}^{-5})$bbNormalization of negative ($n$) and positive ($p$) NPEX components, defined at 1 keV in units of photons keV-1 cm-2 s-1.	${0.5}_{-0.4}^{+1.1}$
$E_{\mathrm{Gaussian}}$ (keV)	19.8(fix)	19.8(fix)
${\sigma }_{\mathrm{Gaussian}}$ (keV)	4.8(fix)	4.8(fix)
$A_{\mathrm{Gaussian}}(\times {10}^{-3})$	${0.5}_{-0.2}^{+0.2}$	${1.0}_{-0.3}^{+0.3}$
$E_{Fe}$ (keV)	6.76(fix)	6.76(fix)
${\sigma }_{Fe}$ (keV)	0.15(fix)	0.15(fix)
$A_{\mathrm{Fe}}(\times {10}^{-4})$ccNormalization of the Gaussian. Defined in units of photons keV-1 cm-2 s-1.	${0.8}_{-0.5}^{+0.5}$	${0.9}_{-0.5}^{+0.5}$
$E_a$ (keV)	—	40.5(fix)
$\sigma$ (keV)	—	${8.5}_{-2.5}^{+3.8}$
$\mathrm{Norm}(\times {10}^{-4})$	—	${2.2}_{-0.7}^{+0.7}$
${\chi }^2$ / d.o.f	395.68/384	366.99/382
$p_{\mathrm{null}}$ddNull hypothesis probability.	3.3$\times$10-1	7.0$\times$10-1

Equations

$$S(E)=\frac{d_{\mathrm{a}}}{\sqrt{2\pi }{\sigma }_{\mathrm{a}}}{e}^{-(E-E_{\mathrm{a}}{)}^2/(2{\sigma }_{\mathrm{a}}^2)},$$

$$I_E(\eta )=N_{\mathrm{p}}\exp \left(-{\left[\frac{\eta -{\bar{\eta }}_{\mathrm{p}}}{\sqrt{2}{\sigma }_{\mathrm{p}}}\right]}^2\right)+N_{\mathrm{f}}\exp \left(-{\left[\frac{\eta -{\bar{\eta }}_{\mathrm{f}}}{\sqrt{2}{\sigma }_{\mathrm{f}}}\right]}^2\right),$$

$${\mathcal{F}}_E(\varphi )=F_1(\varphi ,I_E)+F_2(\varphi ,I_E),$$

$$\Delta =\pi -\arccos \left({\vec{n}}_{AC1}\circ {\vec{n}}_{AC2}\right),$$

Full text

pulsars: individual (4U 1626$-$67) — X-rays: binaries — magnetic fields

\AuthorCallLimit

=20

1 Introduction

The magnetic field strengths of neutron stars can be measured directly by observing cyclotron resonance scattering features (CRSFs, or cyclotron lines) in their hard X-ray spectra. Since the parameters of observed CRSFs are determined by the properties of the accreted plasma (e.g., Mészáros, 1992; Schönherr et al., 2007; Nishimura, 2008; Schwarm et al., 2017a, b), CRSFs in principle provide us with powerful probes of physical processes in strong magnetic fields. This is an active field of research.

The first CRSF was discovered in the X-ray spectrum of Her X-1 by Truemper et al. (1978). To date, CRSFs have been detected from over 20 sources (e.g., Mihara, 1995; Caballero & Wilms, 2012; Tomsick et al., 2015), with observed magnetic field strengths ranging from $10^{12}$ G to $10^{13}$ G. The Nuclear Spectroscopic Telescope Array (NuSTAR; Harrison et al., 2013) is an ideal tool to study cyclotron lines, due to its excellent energy resolution and uninterrupted coverage in the energy band relevant for CRSF discoveries (e.g., Fürst et al., 2014; Tendulkar et al., 2014; Bhalerao et al., 2015; Tsygankov et al., 2016; Bodaghee et al., 2016; Jaisawal & Naik, 2016). On the other hand, modern theoretical models for cyclotron lines are still not fully in agreement with observations. Theory mostly predicts complex line shapes, with simulated CRSFs showing emission wings and asymmetric profiles. Most observations, however, find CRSFs which are well approximated by smooth, symmetric profiles. For example, Her X-1, one of the brightest CRSF sources on the sky, shows no sign of asymmetry or emission wings in its CRSF profile, as shown by Fürst et al. (2013) using NuSTAR and Suzaku data. Legitimate cases of asymmetric CRSFs are rare in the observational literature, with the most notable example probably being Cep X-4, which shows extra absorption in its red wing (Fürst et al., 2015). There is also the case of V 0332+53, which may have an asymmetric profile (Pottschmidt et al., 2005), although NuSTAR data show that the significance of the asymmetric profile depends on the continuum model (Doroshenko et al., 2017).

A notable candidate for a complex CRSF profile is 4U 1626$-$67. This ultracompact X-ray binary (orbital period 42 min; see Middleditch et al., 1981; Chakrabarty, 1998) hosts a 7.7 sec pulsar with a cyclotron line at $\sim$37 keV (Orlandini et al., 1998; Coburn et al., 2002; Iwakiri et al., 2012; Camero-Arranz et al., 2012). The pulsar has undergone two torque reversals in recorded history, in 1990 (Chakrabarty et al., 1997) and 2008 (Camero-Arranz et al., 2010), and is currently spinning up. Suzaku observations bracketing the torque reversal found no changes in the CRSF parameters, despite a factor of $\sim$2.8 increase in X-ray flux and large changes in the soft X-ray spectrum (Camero-Arranz et al., 2012). However, a pulse-phase-resolved study of the 2006 Suzaku observation (during the spin-down state) by Iwakiri et al. (2012) reported the possible detection of an emission-line-like feature at the CRSF energy in the dim phase of the pulse profile. The statistics of the 2010 Suzaku observation were too limited to study the CRSF profile in great detail, despite the higher flux.

An analysis of the NuSTAR observation of 2015 May by D’Aì et al. (2017, hereafter D17) also showed the CRSF to be asymmetric. D17 model the phase-averaged broadband spectrum obtained by NuSTAR and Swift with the bulk+thermal Comptonization continuum model of Becker & Wolff (2007) and an additional component modeled as disk reflection (Ballantyne et al., 2012) with two CRSFs. The profile of the first $37.95\pm 0.15$ keV CRSF is suggested to be complex, and, in contrast to the earlier Suzaku analysis (Camero-Arranz et al., 2012; Iwakiri et al., 2012), a second harmonic at $61.0\pm 1.0$ keV is claimed.

In this paper we re-analyze the NuSTAR and Suzaku data from 4U 1626$-$67 in order to evaluate the significance of the asymmetric line profile of the fundamental CRSF and to perform an analysis of the source behavior using model independent comparisons of the continuum. In addition, to evaluate the geometry of the accretion column, we perform pulse profile modeling using a new relativistic ray tracing code. Comparing the derived geometrical properties and beam patterns with previous theoretical works, we can achieve an understanding of the physical processes in the strong magnetic field. The NuSTAR and Suzaku observations and data reductions are introduced in Section 2.1 and 2.2, respectively. We describe the pulsar’s long-term variability in Section 3.1. In Section 3.2 we model the phase-averaged X-ray spectra using empirical continua and CRSF models. Section 3.3 presents our study of the pulse profile using a new relativistic ray tracing code. Motivated by the pulse profile modeling results, Section 3.4 presents a phase-resolved spectral analysis of the NuSTAR spectrum. In Section 4.1 we discuss the implications and origins of the observed spectral and timing changes between the NuSTAR and Suzaku observations. Section 4.2 investigates the origins of the observed continuum emission. In Section 4.3 we explain the physical interpretation of the pulse profile modeling results comparing with previous theoretical results which take into account the anisotropy of the Thomson scattering cross section in a strong magnetized plasma. Finally, Section 4.4 presents a comprehensive look at the profile of the fundamental CRSF comparing between observations and between pulse phase intervals.

All spectral analysis in this work was performed with XSPEC, v.12.9.0. Unless stated otherwise, all error bars are at the 90% level for one parameter of interest.

2 Observations and Data Reduction

2.1 NuSTAR observation of 4U 1626$-$67

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^{\prime\prime}$ centered on the source.

For the timing analysis, the background light curves were extracted from a circular region with a radius of $100^{\prime\prime}$ 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)111https://github.com/NuSTAR/nuskybgd, using three blank-sky spectra for each telescope extracted from annular regions with radii $300^{\prime\prime}$–$400^{\prime\prime}$, $400^{\prime\prime}$–$500^{\prime\prime}$, and $500^{\prime\prime}$–$600^{\prime\prime}$.

2.2 Suzaku observation of 4U 1626$-$67

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 (CCD) 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 $\mathrm{Gd}_{2}\mathrm{Si}\mathrm{O}_{5}\mathrm{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 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 $4.3^{\prime}$ 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 which 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.

Since the HXD was not capable of imaging, we applied the simulated Non-X-ray Background (NXB) model provided by the Suzaku HXD team (Fukazawa et al., 2009) and a Cosmic X-ray Background (CXB), 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 (SAA) passages and Earth elevation angles result in different exposures for the XIS and the HXD.

2.3 Long-term trend of the X-ray flux below 20 keV

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 $\sim$20% in both the 2–10 and 10–20 keV bands, while the hardness ratio remained constant.

3 Analysis and Results

3.1 Timing analysis

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_{\mathit{Suzaku}}=7.6774(1)$ s and $P_{\mathit{NuSTAR}}=7.67295(1)$ s, respectively. These results are consistent with Fermi/GBM results222https://gammaray.nsstc.nasa.gov/gbm/science/pulsars.html.

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}_{\mathit{Suzaku}}=-2.8\times 10^{-11}\,\mathrm{s}\,\mathrm{s}^{-1}$ and $\dot{P}_{\mathit{NuSTAR}}=-3.3\times 10^{-11}\,\mathrm{s}\,\mathrm{s}^{-1}$, respectively. Using the conservative approach of Takagi et al. (2016), the uncertainty of these $\dot{P}$ values is $0.2\times 10^{-11}\,\mathrm{s}\,\mathrm{s}^{-1}$. 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 years that the period derivative changes are in good agreement with the prediction by Ghosh & Lamb (1979).

3.2 Phase-averaged spectral analysis

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 results333Using the fits for XIS-nominal pointings without XIS2 found in Section 2.2.2 of http://www.astro.isas.jaxa.jp/suzaku/doc/suzakumemo/suzakumemo-2008-06.pdf as a template for the Crab spectrum. We fixed the cross-normalization factor between XIS0 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\pm 0.01$ for XIS0 and $1.13\pm 0.03$ for the HXD-PIN, respectively. Comparing with these results, the cross normalization uncertainty is at the few-percent level.

The CRSF is clearly visible at $\sim$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 $\sim$22.5 keV in the 2010 Suzaku observation to $\sim$21.5 keV in the 2015 NuSTAR observation. A similar spectral difference is also reported by D17, who compared the long-term Swift/BAT 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

[TABLE]

and where $E_{\mathrm{a}}$, $\sigma_{\mathrm{a}}$, and $d_{\mathrm{a}}$ are the energy, width, and line depth. We also apply an iron line modeled by an additive Gaussian, and an additional black body 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}$)444https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl.

The NPEX model succeeded to reproduce the Suzaku data (left of Figure 3a, with $\chi^{2}/\mathrm{d.o.f.}=326.08/317$) but failed to reproduce the NuSTAR data (right of Figure 3a, with $\chi^{2}/\mathrm{d.o.f.}=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 3b 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 3c). A broad Fe K$\alpha$ emission line is also detected at $\sim$6.76 keV with a $\sigma$ of $\sim$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 (FDCO) 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.

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 $\sim$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\pm 1$ keV when using the Becker & Wolff (2007) physical continuum model, compared to their HIGHECUT fits ($67\pm 3$ keV), while the line depth drops from $50^{+17}_{-10}$ to $22\pm 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 241Am 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.6}_{-0.4}$ keV to $5.9\pm 0.1$ keV. While the CRSF is narrower and shallower in the NuSTAR observation compared to Suzaku. The unabsorbed X-ray fluxes in the 3.0–10.0 keV, 30.0–60.0 keV, and 3.0–60.0 keV bands derived with the modNPEX model fitting of the Suzaku data are $3.21^{+0.09}_{-0.12}\times 10^{-10}$, $1.72^{+0.31}_{-0.28}\times 10^{-10}$, and $1.29^{+0.14}_{-0.14}\times 10^{-9}\,\mathrm{erg}\,\mathrm{cm}^{-2}\,\mathrm{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}\,\mathrm{s}^{-1}$ in the 3.0–60.0 keV band. The NuSTAR fluxes in the same energy ranges are $3.44^{+0.04}_{-0.06}\times 10^{-10}$, $1.14^{+0.07}_{-0.08}\times 10^{-10}$, and $1.25^{+0.05}_{-0.04}\times 10^{-9}\,\mathrm{erg}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$, corresponding to an X-ray luminosity of $1.50^{+0.06}_{-0.05}\times 10^{37}\,\mathrm{erg}\,\mathrm{s}^{-1}$ in the 3.0–60.0 keV band. Comparison the observed X-ray fluxes by Suzaku with the NuSTAR results, we found that the 3.0–10.0 keV flux increased by $\sim$5%, whereas the 30.0–60.0 keV flux decreased by $\sim$35%.

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 3c, 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 energies of the GABS models to be the same.

However, since the GABS model only produces an absorption feature, we also tried replacing the second GABS component with an additive Gaussian model (i.e., $[\mbox{modNPEX}+\mbox{Gaussian}]\times\mbox{GABS}$), allowing the Gaussian normalization to be both positive and negative. The best-fit normalization using the model is negative at the $5\sigma$ level, and the energy was consistent with that found using the $2\times\mbox{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\times 10^{-5}$. Therefore, we have possibly detected a complex line profile for the fundamental CRSF in the phase-averaged spectrum.

3.3 Pulse profile modeling

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 eleven 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).

We use the relativistic light bending code of Falkner et al. (2018a, b, see also ) to model the energy-resolved pulse profiles of 4U 1626$-$67 (see Figure 5). This code follows a similar approach to 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.4M_{\odot}$ and radius $R=10$ km. The observed X-rays are emitted by two cylindrical accretion columns AC1 and AC2 of height $h_{\mathrm{AC1,AC2}}$ and radius $r_{\mathrm{AC1,AC2}}$. Allowing for an asymmetric magnetic field, the columns are positioned individually at azimuthal angles $\Phi_{\mathrm{AC1,AC2}}$ and polar angles $\Theta_{\mathrm{AC1,AC2}}$, 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^{\circ}$ 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

[TABLE]

where $\eta$ 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}$, $\sigma$, 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}_{\mathrm{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$, $\sigma$, 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

[TABLE]

where $\phi$ 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}$ (Eq. 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, b).

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$, $\Theta_{\mathrm{AC1,AC2}}$, $\Phi_{\mathrm{AC1,AC2}}$, $h_{\mathrm{AC1,AC2}}$, $r_{\mathrm{AC1,AC2}}$) 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{AC2}}=r_{\mathrm{AC1}}$ and $h_{\mathrm{AC2}}=h_{\mathrm{AC1}}$. Further, the shape of the energy dependent emission profiles of the two columns are tied together, i.e., $\sigma_{\mathrm{p}_{2}}=\sigma_{\mathrm{p}_{1}}$, $\bar{\eta}_{\mathrm{p}_{2}}=\bar{\eta}_{\mathrm{p}_{1}}$ for the pencil beam component, and $\sigma_{\mathrm{f}_{2}}=\sigma_{\mathrm{f}_{1}}$ and $\bar{\eta}_{\mathrm{f}_{2}}=\bar{\eta}_{\mathrm{f}_{1}}=90^{\circ}$ for the fan beam component. We found that there are two local minima of the $\chi^{2}$ landscape of the observer inclination $i$; one solution is $i=9^{\circ}$, the other is $i=27^{\circ}$. Since the $9^{\circ}$ solution shows a simpler pulse profile evolution than the $i=27^{\circ}$ solution and because of its consistency with the physical simulation results taking into account the anisotropy of the scattring 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^{\circ}$ 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 $\chi^{2}$ 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 constraint than the other parameters.

The derived best-fit values of $\Theta_{\mathrm{AC1,AC2}}$, $\Phi_{\mathrm{AC1,AC2}}$ 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 shows that nine of 20 pulsars requires an 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 $\phi=\Phi_{\mathrm{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 $\Theta_{\mathrm{AC1}}-i=3\arcdeg$. Furthermore, the displacement of the two accretion columns compared to the symmetric antipodal case is given

[TABLE]

which is the angular distance of the two unit vectors, $\vec{n}_{\mathrm{AC1},\mathrm{AC2}}$, of the accretion columns positions corresponding to their azimuthal and polar angles. The geometry of our best-fit model yields a displacement of $\Delta=13\fdg 9$, 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 asymmentric 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 an 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 does. This shift explains the widening of the gap between the double peak in the pulse profile. From 3 keV 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 correctness555Otherwise the fit-algorithm might be stuck at a $N_{\mathrm{p}_{2}}=0$ solution, disregarding the second pencil beam solution also for other tested geometries. of our fit we therefore tie the normalization of the two pencil beam components, i.e., $N_{\mathrm{p}_{2}}=N_{\mathrm{p}_{1}}$. Figure 6 shows the energy dependency of the parameter values.

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 since the pencil beam is directed upwards, with only a small offset, $\bar{\eta}_{\mathrm{p}}$, to the magnetic field, which resembles a conical emission pattern. The strong non-pulsed 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}_{\mathrm{p}}$ and beam width $\sigma_{\mathrm{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}_{\mathrm{p_{1}}}$, dominates the evolution of the pulsed fraction, $(\mathcal{F}_{\mathrm{max}}-\mathcal{F}_{\mathrm{min}})/\mathcal{F}_{\mathrm{max}}$, while the width of the fan beam, $\sigma_{\mathrm{f}}$, has only a minor influence (Fig. 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 towards 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, $\eta$, 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, b). As shown in Figure 7, $\eta$ is strongly influenced by the effect of light bending, for the first column $\eta$ varies between $5^{\circ}$ and $22^{\circ}$ in a small band with mean width $\sim 1^{\circ}$, whereas for the second column the mean width of the band is $\sim 11^{\circ}$ between $91^{\circ}$ and $110^{\circ}$.

Compared to 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=9\fdg 5$, 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 $\lesssim 8^{\circ}$ inferred by studies of the orbit of the system that assume that the donor star is a $0.08\,M_{\sun}$ 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_{\sun}$ 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).

3.4 Phase-resolved spectral analysis

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: $\phi=$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 $\eta$ of the first accretion column between $5^{\circ}$ and $16^{\circ}$, while the 0.5–1.0 interval corresponds to angles between $17^{\circ}$ and $22^{\circ}$.

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.

3.4.2 Continuum and cyclotron line

We model the phase-resolved spectra with the modNPEX continuum models, which is the 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 broad Gaussian component for the modNPEX model were fixed at the values from the phase-averaged spectrum (Table 1).

Figure 9 shows the results of the spectral fitting of the spectrum accumulated for phase interval $\phi=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.

We then fitted the spectrum from phase interval $\phi=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 10a), 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 to the GABS CRSF, with $\chi^{2}/\mathrm{d.o.f.}=651.69/427$ (Figure 10b). We then tried another GABS component (i.e., modNPEX$\times$GABS$\times$GABS) and an additive Gaussian component ([modNPEX$+$Gaussian]$\times$GABS). In both cases, the extra residuals around the CRSF were eliminated (Figure 10c,d), and the Gaussian component normalization is negative at the $5\sigma$ level. The probability of this feature arising by chance, determined by simulating 400,000 datasets with simftest in XSPEC, is $1.8\times 10^{-5}$. The best fit parameters using the modNPEX model are summarized in Table 4.

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 ($\phi=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 11a 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 11b), the fitted width of $\sigma=8.5^{+3.8}_{-2.5}$ keV is consistent with the width found in the spin-down state, which showed $\sigma>4.5$ keV, whereas the normalization of $2.2^{+0.8}_{-0.7}\times 10^{-4}$ photons cm*-2* s*-1* is lower than the spin-down observation’s $8.8^{+7.2}_{-2.9}\times 10^{-4}$ ph cm*-2* s*-1*. The chance probability of improvements of adding this emission feature is $2.0\times 10^{-2}$ determined by simulating 20,000 datasets with simftest in XSPEC. Therefore, we concluded that we have only marginally detected the emission line at about 2$\sigma$ 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^{+0.1}_{-1.0}\times 10^{-10}$  $\mathrm{erg}\,\mathrm{cm}^{-2}\,\mathrm{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 $\phi=$0.2–0.35 in Figure 5, and thus this discrepancy is not problematic.

4 Discussion

We have presented an analysis of two observations of 4U 1626$-$67 with Suzaku in 2010 September and NuSTAR in 2015 May. We have performed broad-band spectral analysis, using both phase-averaged and phase-resolved data. In the phase-averaged analysis, we have found a change in the continuum shapes between the observations and confirmed the complex profile of the fundamental CRSF in the NuSTAR data suggested by D17. In the phase-resolved analysis using the NuSTAR data, we have shown the phase dependence of the continuum and CRSF feature, and find further evidence for a distorted CRSF profile during the brighter phase. We have also modeled the energy-resolved pulse profiles using a new relativistic ray tracing code. In this section, we discuss the nature of 4U 1626$-$67, based on our findings.

4.1 Difference of the continuum and pulse period between the two

observations

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\times 10^{-11}$ s s*-1*, while over the NuSTAR observation it was $-3.3\times 10^{-11}$ s s*-1*. 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. Since 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 $\alpha$ 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 to the Suzaku observation. At least in terms of these results, the spectral variation between the 2010 Suzaku and 2015 NuSTAR observations indicate 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.

4.2 Origin of the continuum emission

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 (Ferrigno et al., 2009; Suchy et al., 2008; Vasco et al., 2013; Ballhausen et al., 2017). Since 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 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, since there are three types of seed photons (bremsstrahlung and cyclotron emission created along the column, black-body 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 2nd harmonic absorption feature described earlier. Because of this and since 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.

4.3 Origin of the pulse profile evolution

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 strong 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 May 3–5, 1983. 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 $\tau_{\mathrm{Th}}=80$, Kii et al. (1986) find the following: first, ordinary-mode photons dominate in 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. (1986)’s 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 $\phi=0.45$, where the pencil and fan contributions are comparable, and maximal around phases $\phi=0.3$ and $\phi=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 $\phi=0.7$.

4.4 CRSF profile

We have found that the observed fundamental CRSF in the phase-averaged and the phase $\phi$ = 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 $\phi$ = 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.

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 $\phi$ = 0.0–0.5 spectra, which corresponds to the angle $\eta_{\mathrm{AC1}}$ = 5*∘* – 16*∘* (Figure 7), is well-reproduced by a simple GABS model. Figure 13 shows a comparison of the CRSF shape between phase $\phi$ = 0.5–1.0 and $\phi$ = 0.0–0.5 in terms of data-to-model ratio (with the absorption components excluded). The depth of the line at $\phi$ = 0.0–0.5 is clearly shallower than at $\phi$ = 0.5–1.0. It is noteworthy that the wings of the two profiles resemble each other while the core of the $\phi$ = 0.5–1.0 profile appears to be filled in compared to the $\phi$ = 0.0–0.5 profile. To evaluate the width of the CRSF profiles, we calculate the FWHM (full-width-at-half-maximum) from these ratios. The results are 12.9 keV and 11.4 keV corresponding to $\phi$ = 0.0–0.5 and $\phi$ = 0.5–1.0 , respectively. Therefore, the CRSF is wider and shallower between $\eta_{\mathrm{AC1}}$ = 5*∘* – 16*∘* compared to $\eta_{\mathrm{AC1}}$ = 17*∘* – 21*∘* (see Figure7). The relation between the CRSF shape and the angle $\eta$ is qualitatively consistent with theoretical simulations when assuming a slab-type geometry (Isenberg et al., 1998; Schwarm et al., 2017b). The relation is also qualitatively consistent with the theoretical results assuming a cylindrical geometry illuminated by anisotropic injections and magnetic field gradients (Nishimura, 2015).

The absorption feature disappears in the dim phase spectrum (phase $\phi$=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 $\times 10^{-2}$). The fact indicates that the flux of the emission component decreased to the point of non-detectability 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 4 times lower than during the 2006 Suzaku observation. On the other hand, the 3–60 keV flux of the NuSTAR observation is about 3 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\times 10^{-4}$ using 20,000 datasets with simftest) and the fitted continuum approximating the exponentially cut-off power-law used by Iwakiri et al. (2012) since the contribution of the broad Gaussian and the positive exponential component are comparatively small ($A_{\mathrm{Gaussian}}<0.4\times 10^{-3}$ keV*-1cm-2s-1* and $A_{p}<1.5\times 10^{-5}$ keV*-1cm-2s-1*). According to Section 3.3 the dim phase corresponds to $\eta_{\mathrm{AC1}}\approx 5^{\circ}$, 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.

5 Summary

We have performed a spectral and timing analysis of the accretion powered 7.7 s pulsar 4U 1626$-$67 during its spin-up phase. The results are summarized below:

• •

The Fermi/GBM $\dot{P}$ values during the 2010 Suzaku and 2015 NuSTAR observations are different, implying a $\sim$15% increase in the spin-up rate.

• •

Comparing the phase-averaged 2010 Suzaku and 2015 NuSTAR spectra we found that in addition to the flux increase below 20 keV, the continua differ significantly from each other above 25 keV, with NuSTAR data showing less hard X-ray flux.

• •

Based on the changes in flux, $\dot{P}$, and average spectral shape, we suggest that the accretion rate increased between the Suzaku and NuSTAR observations, associated with decreasing plasma temperature.

• •

Based on the NuSTAR data we confirm earlier results that the pulse profile is strongly energy dependent and changes from being dominated by two narrow peaks in the phase range 0.0–0.5 below 10 keV to being dominated by a single broad peak around pulse phase 0.75 above 20 keV.

• •

The CRSF around 37 keV in the phase-averaged and phase 0.5–1.0 NuSTAR spectra could not be adequately modeled with a single Gaussian (GABS) or pseudo-Lorentzian (CYAB) optical depth profile. A good description was, however, obtained with two GABS components, leading to an asymmetric profile that is shallower towards lower energies.

• •

Possible reasons for an asymmetric CRSF shape are the superposition of line profiles from different locations in the accretion column or the photon spawning effect.

• •

We simultaneously modeled the energy-resolved pulse profiles using a new relativistic ray tracing code to evaluate emission patterns. A combination of pencil- and fan-beam emission with a magnetic field nearly aligned with the line of sight reproduces the data well. In this model the narrow double peak in the pulse profile observed below 10 keV in the first half of the pulse profile is caused by a narrow pencil beam with a small offset to the magnetic field while the flat part of the profile is caused by a fan beam. Towards higher energies the emission geometry of both components evolves (see Figure 5). In particular the offset of the pencil beam regarding the magnetic field enlarges, which causes the double peak to move apart regarding phase increasing its contribution in the second half of the pulse profile.

• •

A comparison of our pulse profile modeling with earlier calculations by Kii et al. (1986) shows that the deduced changes in emission pattern may be due to the energy and polarization dependence of the Thomson scattering cross section in a strong magnetic field.

• •

The CRSF parameters obtained for the two characteristic pulse phase ranges show that the observed CRSF profile depends significantly on the pulse phase. Moreover, our pulse profile modeling leads us to connect the pulse phase with the angle of the emitted photons measured with respect to the magnetic field axis. Connecting these results, we found that the CRSF width decreases and its depth increases with increasing emission angle. This relation is expected from theoretical predictions for slab-type geometries (Isenberg et al., 1998; Schwarm et al., 2017b) as well as in asymmetrically-illuminated cylindrical-geometry models with a magnetic field gradient(Nishimura, 2015).

• •

We also checked for the possible presence of a CRSF signature in emission. Such a feature was tentatively reported for the phase 0.125–0.250 (dim phase) spectrum of the 2006 Suzaku observation, i.e., the spin-down phase, pre-2008-torque-reversal. Such a feature was only marginally detected in the 0.125–0.250 NuSTAR spectrum.

The authors appreciate very much the many constructive comments from the anonymous referee. This work is based on data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the NuSTAR Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). We also thank the Suzaku, Fermi/GBM and MAXI team members for their dedicated support in satellite operations and calibration. This work has been partially funded by the Deutsche Forschungsgemeinschaft under the DFG grant number WI 1860$/$11-1. It uses ISIS functions (ISISscripts) provided by ECAP/Remeis observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/). The figures in this work have been produced with the S-Lang module slxfig. KP acknowledges support by NASA’s NuSTAR Cycle 1 Guest Observer Grant NNX15AV17G. MTW is supported by the Chief of Naval Research. WI is supported by the Special Postdoctoral Researchers Program in RIKEN and JSPS KAKENHI Grant Number 16K17717.