The Modulating Optical Depth of the Photoelectric Absorption Edge with a Pulse Phase in Accretion-powered X-Ray Pulsars

Table 1 summarizes the details of the Suzaku observations of the four selected X-ray pulsars. In all of the observations, the XISs were operated in the normal mode, incorporating a one-fourth window option, which gives a time resolution of 2 s during the observation. In the observation of OAO 1657-415, spaced-row charge injection (SCI) was performed with 2 keV equivalent electrons for XIS 0 and XIS 3 (front-illuminated (FI) CCDs) and 6 keV equivalent electrons for XIS 1 (back-illuminated (BI) CCD), while in the other observations, SCI was conducted with 2 keV equivalent electrons for both FI and BI CCDs. The hard X-ray detector (HXD; Takahashi et al. 2007) was operated in the standard mode, wherein individual events were recorded with a time resolution of 61 μs in all observations.

The total effective exposure times were more than 60 ks in all four observations. In the case of GX 301-2, another Suzaku observation was conducted in 2008 with a total effective exposure time of 10 ks (ObsID=403044010). We did not use these data because of their poor statistics for phase-resolved analysis with a spin period of 700 s.

The archival Suzaku data of the four sources were analyzed using HEASARC software version 6.25. Only XIS 0, 1, and 3 were used and analyzed, since XIS 2 was out of operation due to damage by a micrometeorite in 2006 November. The HXD-PIN and HXD-GSO data were used for broadband spectral analysis. The HXD-PIN data were also used to determine the spin period of the pulsars, since they have better time resolution than XIS. We reprocessed the XIS and HXD data using the FTOOLS task aepipeline (version 2.3.12.25) with the calibration database (versions hxd-20110913, xis-20181010, and xrt-20110630) and the standard screening criteria.

The on-source and background XIS events were extracted from a circular region with a 3' radius and an annulus with an inner radius of 4' and an outer radius of 7' centered at the source position, respectively. For OAO 1657-415 and GX 1+4, a pileup estimation following Yamada et al. (2012) showed that the maximum pileup fraction at the center of the source region is less than 1%, and the effect of the pileup was therefore neglected for these observations. On the other hand, from the on-source events of GX 301-2 and Vela X-1, we excluded the events in central regions with radii of 03 and 09, respectively, to eliminate the pileup effect. This procedure reduces the pileup fraction to below 1%. The redistribution matrix and ancillary response files for each XIS detector were generated using the xisrmfgen and xissimarfgen routines (Ishisaki et al. 2007), respectively. We used tuned (LCFITDT) non-X-ray-background (NXB) models (Fukazawa et al. 2009) to obtain the NXB events for HXD. From the HXD-PIN data, we subtracted the NXB and an expected contribution from the cosmic X-ray background (CXB) with a spectral shape given by Boldt (1987). From the HXD-GSO data, we subtracted only the NXB, as the expected contribution of the CXB is negligible. In our analysis, response files released between 2008 January and 2011 June were used for the HXD-PIN data, selecting a suitable file for each observation, while response and effective area files released in 2010 May were used for the HXD-GSO data.

In the following analyses, all uncertainties quoted are given at the 90% confidence level for the parameter of interest unless stated otherwise.

For the timing analysis, we applied a barycentric correction to the arrival times of individual photons using the aebarycen task FTOOLS (Terada et al. 2008). As for HXD events, dead-time corrections were applied using the FTOOLS task hxddtcor. After the barycentric and dead-time corrections, we conducted a further correction for the effect of the orbital motion in the binary system with the orbital parameters reported by Koh et al. (1997), Kreykenbohm et al. (2008), and Jenke et al. (2012) for GX 301-2, Vela X-1, and OAO 1657-415, respectively. Light curves were extracted from the corrected XIS event data in the range 0.5–12 keV with a 2 s time bin, which is the minimum available time resolution for the XIS operation mode (i.e., one-fourth window option). Data from all of the operated XISs were added together, and a single light curve was obtained. We also created light curves in the range 12–70 keV from HXD-PIN event data with a resolution of 1 s. From the above light curves, we subtracted both the NXB and the CXB.

Figure 2 shows the background-subtracted light curves of the four sources obtained from XIS (upper panel) and HXD-PIN (lower panel), where the time bins are tuned to the pulsar spin period or its harmonic value. A prominent change in the intensity of the light curve of OAO 1657-415 was found during the observation. A low-intensity period was observed during the first ∼1.4 × 10⁵ s, and a high-intensity period followed for ∼6 × 10⁴ s. Since a change of pulse fraction with intensity was reported by Pradhan et al. (2014) for the same Suzaku observation, we divided the data into low- and high-intensity periods at 1.4 × 10⁵ s from the start of the observation and performed separate analyses for the two periods. However, the data in the low-intensity period were not sufficiently statistically sound for a phase-resolved analysis, and we could not extract any meaningful results. In the following description, we show only the results obtained from the last ∼5 × 10⁴ s of the observation, when the source exhibited high intensity. The selected time interval is shown in Figure 2 by the shaded area.

Zoom In Zoom Out Reset image size

To search for pulsations in the light curves of the sources for the phase-resolved analysis, we performed a standard epoch folding analysis (Leahy et al. 1983) of the HXD-PIN data for each observation with the efsearch task of FTOOLS. The barycentric pulsation periods observed in our analysis are summarized in Table 2. Errors estimated according to Leahy (1987) are also given in the table. All of the sources are consistent with the previously reported values derived from the same Suzaku data (Suchy et al. 2012; Maitra & Paul 2013; Pradhan et al. 2014; Yoshida et al. 2017b), within the errors.

Table 2.  Results of Time-averaged Analysis of the Four X-Ray Pulsars from Suzaku Observations

Source	P (s)a	Epochb	Modelc	NH1d	ZFee	NH2f	fpcg	LXh / Di
GX 301-2	685.532(4)	54836.44	FDPL	31 ± 2	1 (fixed)	${33}_{-1}^{+2}$	0.77 ± 0.05	8.5/3.0
Vela X-1	283.449(2)	54634.22	NPEX	${0.77}_{-0.07}^{+0.08}$	12 ± 2j	8.9 ± 0.2	${0.466}_{-0.008}^{+0.007}$	3.1/1.9
OAO 1657-415	36.96(1)	55832.06	NPEX	25.8 ± 0.9	1 (fixed)	${57}_{-5}^{+6}$	${0.59}_{-0.02}^{+0.03}$	20/7.1
GX 1+4	159.944(1)	55471.28	BB+CPL	${18.2}_{-0.4}^{+0.3}$	1.43 ± 0.08	⋯	⋯	8.9/4.3

Notes.

^aBarycentric pulsation period obtained from HXD-PIN data in the range 12.0–70.0 keV. The quoted errors are the statistical 1σ level. ^bEpoch to fold the light curves with individual pulse periods (MJD). ^cApplied continuum model to fit the X-ray broadband spectra. Here BB, CPL, NPEX, and FDPL represent blackbody, exponential cutoff power law, negative and positive power law with an exponential cutoff, and power-law multiplied Fermi–Dirac cutoff, respectively. ^dEquivalent hydrogen column densities (10²² H atoms cm⁻²) estimated from the fully covering photoelectric absorption model. ^eIron abundance derived from the photoelectric absorption model. ^fEquivalent hydrogen column densities (10²² H atoms cm⁻²) estimated from the partially covering photoelectric absorption model. ^gCovering fraction of the partial covering photoelectric absorption model. ^hX-ray luminosity (10³⁶ erg s⁻¹) in the range 0.5–100 keV calculated with a distance of D. Ambiguities of the distance are not included. ⁱDistance to the source from the observer (kpc) from Kaper et al. (2006), Sadakane et al. (1985), Audley et al. (2006), and Hinkle et al. (2006) for GX 301-2, Vela X-1, OAO 1657-415, and GX 1+4, respectively. ^jThis value is thought to be wrong. See text for details.

Download table as:  ASCIITypeset image

Using the estimated pulse periods, the light curves of the XISs were folded by applying the efold task of FTOOLS. Energy-divided pulse profiles obtained from background-subtracted light curves in the 0.5–7.1 and 7.1–12.0 keV ranges are shown in panels (a) and (b) of Figures 3–6, where the data from the three XISs are summed. According to the criteria described in Section 3.3.1, for phase-resolved spectroscopy, the data are divided into pulse phase intervals as indicated by the quoted numbers and overlaying colors in the same panels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows time-averaged and background-subtracted spectra of the four sources, though the spectrum of OAO 1657-415 was obtained from the data of the high-intensity period as described in Section 3.1. The data of Vela X-1 obtained by the three XISs were analyzed in the range 1–10 keV, while for the other cases, we used only 2–10 keV data for FI CCDs. Because they show strong low-energy absorption (more than N_H = 10²³ H atoms cm⁻²), as seen in Figure 1, the detected events below 2 keV are dominated by the "low-energy tail" component (Matsumoto et al. 2006), which is characteristic of the instruments and not well calibrated (Suchy et al. 2012). The XIS data in the 2.2–2.4 keV energy range were ignored for the spectral fitting because of the presence of a gold edge feature, the calibration precision of which is not sufficient compared to the statistical uncertainty. We also omitted the 10–15 keV energy band of the HXD-PIN data because of the calibration uncertainty of the temperature-dependent electrical noise. The cross-normalizations of the XIS detectors were set to be free to cope with the calibration uncertainty of the effective area, whereas the cross-normalizations between XIS and HXD were fixed at 1.16 and 1.18 for "XIS nominal" and "HXD nominal" observations, respectively, following the recommendation in the Suzaku ABC guide.

Using the XSPEC version 12.10 package (Arnaud 1996), these XIS and HXD spectra were simultaneously fit by the model employed by earlier studies to describe the same Suzaku broadband spectra (Suchy et al. 2012; Maitra & Paul 2013; Pradhan et al. 2014; Yoshida et al. 2017b). The spectral continuous components used are summarized in Table 2. A power-law-multiplied Fermi–Dirac cutoff (FDPL; Tanaka 1986) model was applied to GX 301-2, the negative and positive power law with an exponential cutoff (NPEX; Mihara 1995) model was used for Vela X-1 and OAO 1657-415, and an exponential cutoff power law with a blackbody (BB+CPL) model was used for GX 1+4. All of the spectral models were multiplied by the photoelectric absorption TBnew (Wilms et al. 2000) with an improved abundance wilm and cross-section tables vern (Verner et al. 1996). In addition, the model was modified by applying a partial covering intrinsic absorption (e.g., Endo et al. 2000), except for the case of GX 1+4. To represent the cyclotron resonance scattering features in the spectra of GX 301-2 and Vela X-1, the multiplicative components gabs and cyclabs were introduced to the model. In the spectra of Vela X-1 and GX 1+4, we noticed absorption edge–like residuals of data from the model around 7 keV. We fitted this edge-like feature by allowing the iron abundance of the absorber relative to the interstellar medium, Z_Fe, to vary, and values of Z_Fe ∼ 12 in Vela X-1 and ∼1.4 in GX 1+4 were typically obtained. The obtained large abundance of Vela X-1 will be discussed in Section 3.3.1.

The results of the time-averaged spectral fitting are summarized in Table 2. Using the obtained best-fit parameters, we estimated the X-ray luminosities of the four pulsars during the observations (for OAO 1657-415, only during the high-intensity period) listed in Table 2. The distances to the sources from the observer are also given in the table. These values are consistent with those reported by Doroshenko et al. (2011), Suchy et al. (2012), Jaisawal & Naik (2014), and Yoshida et al. (2017b) for the same Suzaku data.

3.3.1. Spectroscopy

We performed phase-resolved spectroscopic analyses of the four sources, focusing on the emission lines and absorption edge of iron. To study the variation of the spectral shape according to the pulse phase, we divided the observation data into several phase bins, as indicated in Figures 3–6 with different overlaid colors. The bins were chosen to pick up various prominent features, e.g., the deep minimum and maximum peaks, as well as to retain photon statistics sufficient to constrain the individual spectral parameters. The hardness ratios between the two pulse profiles plotted in panel (c) of Figures 3–6 were examined, and the boundaries were found to be correct for picking up a stable hardness ratio region, except for continuously changing phases. In the case of GX 1+4, the phase intervals 1 and 3 (marked in purple in Figure 6) were combined to improve the photon statistics.

We fitted the broadband phase-resolved spectra with the same models as those used in the time-averaged spectral fitting listed in Table 2. When the phase-resolved spectra of GX 301-2 were fitted using the FDPL model, we fixed the cutoff energies at 29 keV, which is the best-fit value of the time-averaged spectrum, whereas the other parameters were left free. In the fitting of Vela X-1 and GX 1+4, Z_Fe was fixed at the value obtained from the time-averaged fitting. We found that the adopted model fit all of the phase-resolved spectra well.

The best-fit parameters of the column density of the photoelectric absorption (of the fully covering component for GX 301-2, Vela X-1, and OAO 1657-415) and the photon index (of the CPL model, the negative power-law component in the NPEX model, and the FDPL model) derived from the above fitting are plotted as a function of the pulse phase in panel (d) of Figures 3–6. The figures show strong positive correlations between the column density and the photon index suggesting that the derived values influence each other. For example, if the spectral shape is steep in the high-energy band, the derived column density value becomes high because the extrapolation to the low-energy band of the steep spectrum requires a large photoelectric absorption, without any evidence to justify the extrapolation. Therefore, the absorption column densities are not physically meaningful values, since we cannot apply any physically well-justified model in the low-energy region. If the column density is incorrect, the K_α line flux and the K-edge optical depth of the iron ions are also affected, and we may derive wrong values. We consider that the iron abundance of 12 ± 2 of Vela X-1 obtained from the time-averaged spectral fit, shown in Table 2, is wrong.

For a robust analysis of the variation of the emission lines and absorption edge with the pulse phase, avoiding the influence from the insufficient modeling of the broadband spectrum, we restricted the energy range to 5.0–7.9 keV for the spectral fitting. In this restricted region, the continuum can be fitted by a power-law model multiplied by an absorption edge; hence, we can deduce the edge parameter independently from the low-energy absorption. An approximate model, edge, is employed to reproduce the photoelectric absorption edge, which is given as follows:

$$\begin{eqnarray}M(E)=\left\{\begin{array}{ll}1 & (E\lt {E}_{{\rm{c}}})\\ \exp \left(-{\tau }_{\mathrm{edge}}{\left(\displaystyle \frac{E}{{E}_{{\rm{c}}}}\right)}^{-3}\right) & (E\geqslant {E}_{{\rm{c}}}),\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where E_c is the threshold energy and τ_edge is the absorption optical depth at the threshold energy. The applied model consists of powerlaw × edge + Gaussian + Gaussian + Gaussian, in which three Gaussian models express the iron K_α, iron K_β, and nickel K_α emission lines. Since the iron K_β line contaminates the absorption edge feature, it is difficult to independently determine the parameters corresponding to these structures. Since the energy of the observed iron K_α emission line is around 6.4 keV, the ionization state of the iron should be low (Fe_XVII at most). In the above ionization state, the difference of the ratio between the energies of the iron K_α and iron K_β lines is at most 1%, but that of the energy of the iron edge to that of the iron K_α line is at most 7% (Kallman et al. 2004; Yamaguchi et al. 2014). Therefore, we fixed the ratio of the energy of the iron K_β line to that of the iron K_α line at 1.103 for the neutral case, while the energy of the iron K-edge was allowed to vary in the fitting. In the case of OAO 1657-415, as the center energy of the nickel K_α line was poorly constrained due to low effective exposure in each phase bin, the ratio of its energy was fixed to that of the iron K_α line at 1.168 for the neutral case. The widths of the emission lines could not be resolved in the XIS data, and they were fixed to be sufficiently narrow. Note that for our broadband spectroscopy, we used only the FI CCDs for GX 301-2, OAO 1657-415, and GX 1+4 because of a discrepancy between the FI and BI CCDs in the low-energy region (below 3 keV). However, in the above restricted narrow energy range, the discrepancy between the FI and BI CCDs is negligibly small, and we therefore used data from all of the CCDs.

Parameters obtained from the spectral fitting are summarized in Table 3 and plotted as a function of pulse phase in Figures 3–6. We find that the optical depths of the iron K-edge, τ_edge, exhibit modulations with pulse phase, as seen in panel (e) of each figure. For all four pulsars, τ_edge increases at the phase when the X-ray continuum flux becomes dim. The energies of the iron K_α line and iron K-edge remain constant within the errors, as shown in panels (f) and (g). The ratios of the energy of the iron K-edge to that of the iron K_α line are calculated and plotted against the pulse phase in panel (h), and the values are constant within the errors. Note that for OAO 1657-415 and GX 1+4, the energy ratio shows marginal modulation due to the change of the energy of the edge, though with poor statistical significance, which peaks synchronously with the maximum τ_edge. The expected energy ratios for Fe_I–Fe_VII are also shown in the same panels by gray dashed lines. The averaged values of the energy ratios of the four pulsars are in the range 1.117–1.122. These values show that the ionization state of the iron is Fe_II–VII. We confirmed that the assumed ratio of the energy of the iron K_β line to that of the iron K_α line of 1.103 in the fitting is consistent with this ionization state. Table 3 lists the obtained intensities of the iron K_α and K_β. The ratios between these intensities are roughly consistent with the expected values in the above ionization state of 0.12–0.14 (Yamaguchi et al. 2014), except for the results of GX 301-2 and interval 3 of OAO 1657-415. In these cases, the assumed energy ratio between the iron K_α and K_β lines might be smaller than the actual value; thus, the results for the edge parameters of GX 301-2 and OAO 1657-415 might be marginal. However, in the case of Vela X-1 and GX 1+4, the edge parameters are not affected by the contamination of the iron K_β line.

Table 3.  Results of Phase-resolved Spectral Fitting

Intervala	${E}_{{\mathrm{FeK}}_{\alpha }}$ b	${E}_{{\mathrm{FeK}}_{\mathrm{edge}}}$ b	τedgec	NHd	${I}_{{\mathrm{FeK}}_{\alpha }}$ e	${I}_{{\mathrm{FeK}}_{\beta }}$ e	${\chi }_{\nu }^{2}$ (d.o.f.)
GX 301-2
1	6.392 ± 0.004	${7.152}_{-0.01}^{+0.007}$	${0.68}_{-0.04}^{+0.02}$	6.4	${2.04}_{-0.07}^{+0.09}$	$\lt 0.08$	0.98(1663)
2	6.392 ± 0.002	${7.142}_{-0.01}^{+0.007}$	0.72 ± 0.03	6.8	1.97 ± 0.05	<0.08	1.09(2126)
3	${6.392}_{-0.005}^{+0.004}$	${7.145}_{-0.007}^{+0.006}$	0.65 ± 0.02	6.2	2.18 ± 0.08	<0.06	1.14(2241)
4	${6.392}_{-0.003}^{+0.002}$	${7.147}_{-0.007}^{+0.006}$	0.67 ± 0.02	6.3	${2.34}_{-0.08}^{+0.07}$	<0.02	1.15(2274)
5	${6.392}_{-0.004}^{+0.005}$	${7.151}_{-0.01}^{+0.008}$	${0.77}_{-0.04}^{+0.03}$	7.3	${2.17}_{-0.07}^{+0.08}$	<0.09	1.09(1355)
6	6.394 ± 0.004	${7.155}_{-0.009}^{+0.008}$	0.70 ± 0.03	6.6	2.07 ± 0.08	<0.03	1.01(1769)
Vela X-1
1	6.403 ± 0.003	7.17 ± 0.03	0.20 ± 0.02	1.9	1.88 ± 0.07	${0.09}_{-0.08}^{+0.09}$	1.00(2329)
2	6.400 ± 0.004	7.21 ± 0.03	0.17 ± 0.02	1.6	1.71 ± 0.09	<0.2	0.99(2373)
3	${6.405}_{-0.007}^{+0.005}$	${7.17}_{-0.05}^{+0.04}$	0.12 ± 0.02	1.1	1.7 ± 0.1	<0.2	1.02(2323)
4	6.402 ± 0.006	${7.19}_{-0.05}^{+0.04}$	0.11 ± 0.02	1.0	1.8 ± 0.1	0.2 ± 0.1	0.99(2312)
5	${6.402}_{-0.002}^{+0.004}$	${7.17}_{-0.03}^{+0.02}$	0.17 ± 0.02	1.6	1.79 ± 0.06	${0.07}_{-0.07}^{+0.08}$	1.03(2377)
6	${6.400}_{-0.004}^{+0.002}$	7.16 ± 0.04	0.16 ± 0.02	1.5	1.82 ± 0.06	${0.1}_{-0.09}^{+0.1}$	1.04(2377)
OAO 1657-415
1	6.410 ± 0.006	7.22 ± 0.02	${0.75}_{-0.06}^{+0.09}$	7.1	1.17 ± 0.06	0.07 ± 0.05	1.15(458)
2	${6.412}_{-0.006}^{+0.004}$	${7.21}_{-0.03}^{+0.02}$	0.63 ± 0.05	6.0	1.23 ± 0.06	<0.13	1.05(583)
3	${6.410}_{-0.004}^{+0.006}$	${7.20}_{-0.02}^{+0.01}$	0.61 ± 0.04	5.8	${1.23}_{-0.03}^{+0.06}$	<0.05	0.99(853)
4	${6.408}_{-0.004}^{+0.005}$	7.19 ± 0.02	${0.64}_{-0.04}^{+0.05}$	6.1	1.21 ± 0.06	${0.09}_{-0.08}^{+0.07}$	1.02(772)
5	${6.410}_{-0.004}^{+0.006}$	${7.18}_{-0.03}^{+0.02}$	0.57 ± 0.05	5.4	${1.16}_{-0.07}^{+0.06}$	<0.1	0.99(619)
6	${6.409}_{-0.007}^{+0.005}$	${7.19}_{-0.03}^{+0.02}$	0.69 ± 0.06	6.5	1.15 ± 0.06	<0.12	1.02(616)
GX 1+4
1 and 3	${6.427}_{-0.004}^{+0.005}$	7.21 ± 0.02	0.43 ± 0.04	4.1	1.11 ± 0.05	<0.1	0.98(1351)
2	${6.422}_{-0.003}^{+0.004}$	${7.23}_{-0.01}^{+0.02}$	0.52 ± 0.06	4.9	${1.12}_{-0.04}^{+0.05}$	${0.06}_{-0.05}^{+0.04}$	0.93(889)
4	6.425 ± 0.005	${7.20}_{-0.03}^{+0.02}$	${0.37}_{-0.04}^{+0.03}$	3.5	${1.02}_{-0.04}^{+0.05}$	${0.10}_{-0.06}^{+0.05}$	0.94(1900)
5	${6.424}_{-0.002}^{+0.004}$	${7.21}_{-0.03}^{+0.02}$	0.37 ± 0.03	3.5	1.02 ± 0.04	$\lt 0.07$	0.97(2172)
6	${6.422}_{-0.003}^{+0.004}$	7.20 ± 0.02	0.35 ± 0.03	3.3	1.07 ± 0.04	$\lt 0.06$	1.01(2212)
7	6.426 ± 0.004	${7.19}_{-0.03}^{+0.02}$	0.36 ± 0.03	3.4	1.15 ± 0.05	${0.11}_{-0.07}^{+0.08}$	0.90(2027)
8	${6.423}_{-0.002}^{+0.005}$	7.19 ± 0.03	0.37 ± 0.04	3.5	${1.16}_{-0.05}^{+0.06}$	0.2 ± 0.1	0.93(1570)

Notes.

^aPulse phase intervals corresponding to the quoted numbers and different overlaying colors in Figures 3–6. ^bEnergies of the K_α emission line and absorption K-edge of iron in units of keV. ^cOptical depth at the absorption edge of iron. ^dEquivalent hydrogen column density (10²³ H atoms cm⁻²) calculated from the optical depth τ_edge assuming the abundance of the solar system (Asplund et al. 2009). ^eFluxes of the iron K_α and K_β emission lines in units of 10⁻³ photons s⁻¹ cm⁻².

Download table as:  ASCIITypeset image

Using the obtained optical depth, an equivalent hydrogen column density for an assumed solar abundance of ${N}_{{\rm{H}}}=\tfrac{{\tau }_{\mathrm{edge}}}{{Z}_{\mathrm{Fe}}{\sigma }_{\mathrm{Fe}}}$ is derived, where Z_Fe and σ_Fe are the iron abundance in the solar system (=3.2 × 10⁻⁵; Asplund et al. 2009) and the photoelectric absorption cross section of Fe_II at the energy of the K-edge (=3.3 × 10⁻²⁰ cm²; Verner & Yakovlev 1995). The estimated equivalent hydrogen column densities for each phase bin are given in Table 3.

3.3.2. Investigation of Statistical Significance of Modulation

As noted in the above section, modulations of the optical depth of the iron K-edge with pulse phase can be seen in all four pulsars. We performed χ² tests to evaluate the statistical significance, with a null hypothesis that there is no variation of the optical depth over the pulse phase. The χ² values obtained from the constant fittings to the optical depth in the four pulsars are summarized in Table 4 with the corresponding null hypothesis probabilities. These obtained values clearly indicate that the pulse phase modulations in the optical depth are statistically significant. To further inspect the significance of the variation in the optical depth, we directly compared the spectra extracted from different phase bins by calculating their spectral ratio. For clarity, two phase-resolved spectra were chosen so as to maximize the optical depth difference. If the optical depths are different and the energies of the edge are the same between the two phases, a stepwise feature is expected around the energy of the edge in the spectral ratio. The calculated spectral ratio is shown in panel (a) of Figures 7–10. In the spectral ratios, edge features can be seen around 7.1 keV. In addition, excesses are exhibited at 6.4 and 7.5 keV, resulting from variations between the two phases of the equivalent widths of the iron K_α line and nickel K_α line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Summary of the Statistical Test to Evaluate the Modulation in the Optical Depth of the Iron K-edge with Pulse Phase

Constant Fitting		Spectral Ratio Fitting					Run Test
χ2(d.o.f.)a	Probabilitya	Δ	${\chi }_{1}^{2}({\nu }_{1})$ b	${\chi }_{2}^{2}({\nu }_{2})$ b	Fc	Probabilityc	K/n+/n−d	Probabilitye
GX 301-2
44.16(5)	2.14 × 10−8	0.11 ± 0.03	91.55(52)	73.33(51)	1.22	2.35 × 10−1	9/15/11	2.71 × 10−2
Vela X-1
33.33(5)	3.24 × 10−6	0.10 ± 0.02	87.22(52)	55.54(51)	1.54	6.25 × 10−2	5/10/16	2.15 × 10−4
OAO 1657-415
18.10(5)	2.83 × 10−3	0.21 ± 0.04	60.10(52)	41.77(51)	1.41	1.10 × 10−1	8/17/9	1.71 × 10−2
GX 1+4
30.17(7)	8.84 × 10−5	0.20 ± 0.03	88.35(52)	54.42(51)	1.59	4.93 × 10−2	8/12/14	8.52 × 10−3

Notes.

^aThe χ² value with degree of freedom (quoted value) and corresponding probability for the null hypothesis that there is no variation. ^bTwo χ² values corresponding to the results of fittings using models without or with the edge feature component. The quoted values are degrees of freedom. ^cThe F statistical value and corresponding probability of chance improvement of χ² with an F-test. ^dNumber of runs and data points above and below zero. ^eNull hypothesis probability of the randomness in the residuals of the spectral ratio fitting.

Download table as:  ASCIITypeset image

Next, we fitted the spectral ratio in the 5.0–8.4 keV range with the following empirical model:

$$\begin{eqnarray}&&R(E)=\left(A\times E+B\right)\times M(E)+{G}_{\mathrm{Fe}}(E)+{G}_{\mathrm{Ni}}(E),\end{eqnarray} \tag{ 2 }$$

where A is the slope, B is the intercept, M(E) is the edge component given by Equation (1), and G_Fe(E) and G_Ni(E) are Gaussian functions expressing the emission line features at 6.4 and 7.5 keV. In this case, the τ_edge in M(E) roughly corresponds to the difference in the optical depth at the iron K-edge between the two spectra, denoted by Δ. In our fitting, the energy of the edge feature was fixed to the value obtained from the phase-resolved spectral fitting at the phase when the optical depth is maximum. The feature of the nickel emission line around 7.5 keV is not clearly separated from the edge feature in the spectral ratio. Therefore, the ratio of the center energy of G_Ni(E) to that of G_Fe(E) is fixed at 1.168 for the neutral case, and the width of G_Ni(E) is assumed to be equal to that of G_Fe(E). The other parameters were allowed to vary. In order to simulate smearing due to the detector response of XIS, the model was convolved with a single Gaussian function with a width representing the energy resolution at the observation period.

The best-fit models are plotted in panel (a) of each figure, and the residuals from the best-fit model are plotted in panel (b). Panel (c) shows residuals from the best-fit model without introducing the edge component, showing stepwise features around 7 keV. The depths of the edge feature Δ derived from the fitting are listed in Table 4. It is confirmed that the values of Δ are consistent with the differences between the optical depths of the two phase-resolved spectra, within the errors, for all sources (see Table 3). The resultant χ² values are also summarized in Table 4.

To evaluate the statistical significance of the χ² improvement due to the addition of the edge component, we performed an F-test routine described in Press et al. (2007). The F statistical value is defined as

$$\begin{eqnarray}&&F=\displaystyle \frac{{\chi }_{1}^{2}/{\nu }_{1}}{{\chi }_{2}^{2}/{\nu }_{2}},\end{eqnarray} \tag{ 3 }$$

where ${\chi }_{1}^{2}$ and ${\chi }_{2}^{2}$ are χ² values and ν₁ and ν₂ are degrees of freedom corresponding to the results of fittings using models 1 and 2 (in this case, model 1 is a linear function along with the two line components, and model 2 is a linear function multiplied by the edge component with the two line components), respectively. The derived F statistical values are listed in Table 4. The corresponding chance probabilities in the random case of an improvement of χ², given in Table 4, are not small, and we can only marginally conclude the differences of the optical depths by the F-test.

However, the residuals shown in panel (c) of the above figures clearly show edge features around 7 keV. The edge feature is like a stepwise function and might be sensitively tested by a run test (also called the Wald–Wolfowitz test). Therefore, we applied the run test to the residuals obtained from the spectral ratio fitting using the linear function along with the two line components. We encountered a technical problem in this procedure, however, in that we obtained the wrong normalization of G_Ni(E) in the fitting without the edge component. When we fit the data without the edge component, the best-fit continuum component (linear component) became more intense than the data in the energy range above the edge energy and faint below it. Even if a hump due to the nickel line exists in the data, the model simulating the nickel line can only reproduce the excess above the continuum model, which should already be stronger than the data. Therefore, the obtained best-fit normalization of the nickel line by the model without the edge component is smaller than the actual value in the data. As a result, the G_Ni(E) feature remains in the residuals and prevents a correct evaluation of the edge feature with the run test. Therefore, for the run test, we use the residuals produced by the fitting without the edge component whose normalization of G_Ni(E) is fixed at the value obtained from the fit using the edge component. The residuals are plotted in panel (d) of each figure. The run test was conducted to evaluate the null hypothesis of the randomness in the residuals in the range 6.5–8.0 keV. The number of runs and data points below and above zero is summarized in Table 4, and the derived null hypothesis probabilities from the above run test are also given. In the cases of Vela X-1 and GX 1+4, the computed probability of random sampling is <0.01. These probabilities suggest a significant detection of the edge-like features in the spectral ratio; hence, the optical depth is indeed different between the two phase bins. For GX 301-2 and OAO 1657-415, similar edge-like features can be seen in the residuals, although they have poor statistical significance from the run test.

Below, we discuss the following four possible interpretations to explain the modulating optical depth of the absorption edge with the pulse phase: (1) a variation in the physical state of the absorbing matter, such as the ionization state, according to the spin phase; (2) the existence of a modulating component reflected by some neutral matter; (3) a movement of the emission region; and (4) a change of the absorbing matter location.

4.1.1. Physical State Variation of Absorbing Matter

4.1.2. Modulating Reflected X-Rays

4.1.3. Movement of the Emission Region

4.1.4. Change of Absorbing Matter Location

Alternatively, if the absorbing matter corotates with the neutron star spin, variation of the optical depth with the pulse phase can occur. The corotation with the neutron star spin implies that the absorbing matter should be confined by the magnetic field of the neutron star; hence, it should be located at or within the Alfvén radius R_A. The outer accretion disk and wind from companion stars cannot be the cause of the pulse phase modulation of the optical depth, because these locations cannot have variations with the spin period. The Alfvén radius is given as

$$\begin{eqnarray}&&{R}_{{\rm{A}}}=3.7\times {10}^{8}{\left({M}_{* }/{M}_{\odot }\right)}^{1/7}{R}_{6}^{10/7}{B}_{12}^{4/7}{L}_{37}^{-2/7}\ \mathrm{cm},\end{eqnarray} \tag{ 4 }$$

where M_*, R₆, B₁₂, and L₃₇ are the mass and radius (denoted as R_*) in units of 10⁶ cm, surface magnetic field strength (denoted as B_surf) in units of 10¹² G of the neutron star, and luminosity (denoted as L_X) in units of 10³⁷ erg s⁻¹, respectively (Frank et al. 2002). Substituting a neutron star radius of R_* = 10 km and a mass of M_* = 1.4 M_⊙ in Equation (4), the Alfvén radii for Vela X-1 and GX 1+4 can be calculated using the estimated X-ray luminosity in the Suzaku observations. Here the magnetic field strength at the neutron star surface of Vela X-1 is assumed to be 2.2 × 10¹² G, derived from the cyclotron resonance scattering feature in the Suzaku X-ray spectrum (Doroshenko et al. 2011). However, we tentatively assume that of GX 1+4 to be 10¹³ G, since no cyclotron resonance scattering features have been reported in its spectra (Yoshida et al. 2017b).

The absorbing matter responsible for the absorption edge should be located within the Alfvén radius. This fact puts restrictions on the distance between the X-ray source and the matter, r, as

$$\begin{eqnarray}&&r\leqslant 8.1\times {10}^{8}\,\mathrm{cm}\end{eqnarray} \tag{ 5 }$$

for Vela X-1 (L_X = 3.1 × 10³⁶ erg s⁻¹ and B_surf = 2.2 × 10¹² G) and

$$\begin{eqnarray}&&r\leqslant 1.4\times {10}^{9}\,\mathrm{cm}\end{eqnarray} \tag{ 6 }$$

for GX 1+4 (L_X = 8.9 × 10³⁶ erg s⁻¹ and B_surf = 1.0 × 10¹³ G). The absorbing matter should also satisfy the constraint due to the ionization parameter based on the observed ionization state of iron. Specifically, the restrictions are written as

$$\begin{eqnarray}&&{nr}^{2}\gt 6.9\times {10}^{34}\,{{\rm{c}}{\rm{m}}}^{-1}\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{{nr}}^{2}\gt 2.0\times {10}^{35}\,{\mathrm{cm}}^{-1}\end{eqnarray} \tag{ 8 }$$

for Vela X-1 and GX 1+4, respectively, where n is the number density of the particles of the matter. The acceptable regions on the r–n planes satisfying the requirement on the absorbing matter given by Equations (5)–(8) are shown in light gray in panels (a) and (b) of Figure 11 for Vela X-1 and GX 1+4, respectively. For the observed amount of change over the pulse phase in the hydrogen column density of Δ N_H = 10²³ cm⁻², the geometric thickness of the absorbing matter along the line of sight, δ, can be calculated for a given particle density and is indicated by the vertical axes on the right of the panels of Figure 11. If the absorbing matter is at the Alfvén radius and has a particle density of n = 10¹⁸ cm⁻³, the absorbing matter contains almost neutral iron and forms a structure whose geometric thickness along the line of sight is 10⁵ cm. Given these facts, the observed variation in the amount of absorption in the line of sight due to the absorbing matter corotating with the neutron star spin is a plausible explanation for the modulating optical depth of the iron absorption edge with the pulse phase.

Zoom In Zoom Out Reset image size

The condition described in Section 4.1.4 suggests that the absorbing matter is the matter accreting onto the neutron star along its magnetic field lines. We can therefore discuss the geometry of the accretion flow and determine the relation between the mass accretion rate $\dot{M}$, r, n(r), and the falling velocity v(r). We thus obtain

$$\begin{eqnarray}\begin{array}{rcl}\dot{M} & = & {\mu }_{\mathrm{gas}}{m}_{{\rm{p}}}n(r)s(r)v(r)\\ & = & \sqrt{2{{GM}}_{* }/r}{\mu }_{\mathrm{gas}}{m}_{{\rm{p}}}n(r)s(r),\end{array}\end{eqnarray} \tag{ 9 }$$

where s(r) is the cross-sectional area of the accretion flow, m_p is the proton mass, and μ_gas is the mean atomic mass of the cosmic material, which is 1.3 for solar abundances (Cox 2000). We assumed the velocity falling onto the neutron star to be the freefall velocity given by $v(r)=\sqrt{2{{GM}}_{* }/r}$. By assuming that all of the kinetic energy liberated is converted to radiation energy at the neutron star surface, namely, ${L}_{{\rm{X}}}={{GM}}_{* }\dot{M}/{R}_{* }$, the density in the accretion flow is given by

$$\begin{eqnarray}\begin{array}{rcl}n(r) & = & 2.5\times {10}^{27}{r}^{1/2}s{\left(r\right)}^{-1}\\ & & \times \,{\left({M}_{* }/{M}_{\odot }\right)}^{-3/2}{R}_{6}{L}_{37}\ {\mathrm{cm}}^{-3}.\end{array}\end{eqnarray} \tag{ 10 }$$

The cross-sectional area of the accretion flow is given by the product of the width and thickness of the accretion flow, w(r) and d(r), and is assumed to be proportional to r^γ. A schematic picture of the assumed accretion flow is shown in Figure 12. We assume that the width is w(r) = ζr, where ζ is a constant azimuthal angle of the accretion flow, and that the thickness can be written as $d(r)={d}_{{\rm{A}}}{\left(r/{R}_{{\rm{A}}}\right)}^{\gamma -1}$, where d_A is the thickness at R_A. Then, we obtain

$$\begin{eqnarray}&&s(r)=\zeta {d}_{{\rm{A}}}{{R}_{{\rm{A}}}}^{1-\gamma }{r}^{\gamma }.\end{eqnarray} \tag{ 11 }$$

In the case of GX 1+4, the azimuthal angle of the accretion flow can be estimated from the pulse phase duration during the deepest edge phase, yielding ζ = 0.4π (see Figure 6). We also assume the same condition for Vela X-1, although we have no basis for this assumption. If we consider a simple conic geometry for the accretion flow, the index of γ is 2. On the other hand, Ghosh & Lamb (1979) reported that the cross-sectional area of the accretion flow s(r) can be assumed to be proportional to the cube of r, namely, γ = 3. Since this assumption of γ = 3 can be applied only near the neutron star and cannot be extrapolated to the Alfvén radius, we consider that the actual γ value should be between 2 and 3. Thus, we calculated cases of both γ = 2 and 3. Substituting Equation (11) into Equation (10) with γ = 2 and 3, we calculated the number density in the accretion flow. Using the obtained parameters for Vela X-1 and GX 1+4, i.e., the X-ray luminosities, magnetic strengths at the neutron star surface, and ζ = 0.4π, the calculated number densities in the accretion flow are plotted by red dashed lines (γ = 2) and blue dotted lines (γ = 3) in Figure 11 for some typical thicknesses for the accretion flow at the Alfvén radius, d_A = 10⁴, 10⁵, and 10⁶ cm. For both Vela X-1 and GX 1+4, as seen in Figure 11, the accreting matter at any radius is acceptable if it has a thickness of d_A < 10⁵ cm and its cross-sectional area is proportional to r³. On the other hand, for γ = 2, only the outer region is acceptable for accreting matter of thickness d_A ∼ 10^4–5 cm for both sources. It should be noted that the thickness of the modeled accretion flow at the Alfvén radius, d_A, is roughly consistent, within a factor of 3–8, with δ, indicated on the vertical axes on the right of the panels in Figure 11.

Zoom In Zoom Out Reset image size

As mentioned in Section 3.3.1, the deepest edge phase corresponds to the dimming interval of the X-ray continuum, such as a dip interval. The dip in the pulse profile is interpreted as being due to the eclipse of the X-ray-emitting region by the accretion column of the pulsar, at least for GX 1+4 (Dotani et al. 1989; Galloway et al. 2000, 2001; Giles et al. 2000). Galloway et al. (2001) reported that the absorption column density, derived from the low-energy cutoff, increased at the dip. This is in good agreement with the results described here. The increasing amount of absorption could be attributed to a crossing of the line of sight by the accretion column or curtain (Miller 1996). In our geometry of the accretion flow with Equation (11), its thickness is sufficiently smaller than its width if the thickness is the estimated value of 10^4–6 cm at the Alfvén radius. Therefore, for the configuration of the accretion flow, a geometry that is flattened along the azimuthal direction like an accretion curtain (Miller 1996; Galloway et al. 2001; the latter shows the schematic figure of accretion geometry) may be more plausible than a whole cone geometry.

The observed modulating optical depth of the iron K-edge with pulse period can be explained by the absorption by the accreting matter along the magnetic field lines within the Alfvén radius, which corotates with the neutron star. That is to say, the shallowest edge phase can be interpreted as the phase when the accreting matter is out of the line of sight. However, we note that even at the shallowest edge phase, a large amount of absorption, corresponding to N_H = 10²³ cm⁻², was observed (see Table 3). We cannot determine whether this large absorption is caused by the accreting matter within the Alfvén radius or by the circumstellar matter located outside the Alfvén radius.

Several authors have discussed the possibility that a reprocessing site contributes to the 6.4 keV iron fluorescence line emission. For Vela X-1, the presence of a line emission site within 5 × 10¹¹ cm of the neutron star was proposed as an explanation for a drop of the iron line flux of a factor of 20 during the eclipse observed by Ginga (Ohashi et al. 1984). For GX 1+4, it was reported that the peak phase of the modulating line flux is delayed compared to that of the pulse profile of the continuum by ∼30 s (Yoshida et al. 2017b). A possible explanation for this delay was proposed by Yoshida et al. (2017a), where the iron line is a fluorescence of the matter exposed to the X-ray continuum from the neutron star at a distance of ∼10¹² cm from the X-ray source. Therefore, it is natural that a certain amount of matter exists outside the Alfvén radius. If so, the matter located in this region should be responsible for not only the fluorescent but also the photoelectric absorption. The absorption contribution by the matter in this region is not affected by the spin of neutron stars and is a possible candidate for the absorption observed at the shallowest edge phase.

In the above section, we argued that the absorbing matter may be distributed in two regions separated by the Alfvén radius. The matter within the Alfvén radius has an asymmetric structure and is the origin of the modulating absorption, while the matter outside the Alfvén radius is responsible for the phase-independent absorption. In our analysis described in Section 3.3, we assumed that the absorbing matter is one component, and that there is no variety in the ionization state of the absorbing matter along the line of sight. We then determined the ionization state of iron from the ratio between the energies of the iron K_α emission line and the iron absorption K-edge. However, the ionization states of the absorbing matter in the two regions does not need to be the same. The matter causing the modulation in the absorption and the other matter causing constant absorption should be handled separately in the analysis. Therefore, we conducted an additional spectral analysis dealing with the absorption by the matter in the two regions separately.

We fitted the phase-resolved spectra of Vela X-1 and GX 1+4 at the deepest edge phase, restricting the energy range to 5.0–7.9 keV, similar to that in Section 3.3. For the spectral fitting in the restricted energy range, a model consisting of powerlaw × edge₁ × edge₂ + Gaussian + Gaussian + Gaussian is applied to express the absorption by iron in the two regions. In this model, edge₁ and edge₂ represent the absorption edges of the matter outside and within the Alfvén radius, respectively. In the fitting, the ratio of the energy of the iron K_β line to that of the iron K_α line was fixed to 1.103 for the neutral case (Yamaguchi et al. 2014). The center energy of the nickel K_α line was also fixed to a value derived from the phase-resolved fitting at the deepest edge phase. The energy and optical depth of the edge₁ component were fixed to the values obtained from the fitting of the shallowest edge phase. The energy of the edge₂ component was scanned by changing the ratio to the fixed edge energy from 1.00 to 1.05 with steps of 0.001, while the optical depth of the edge₂ component was allowed to be free. The resultant χ² values as a function of the energy of edge₂ are plotted in the upper panels of Figures 13(a) and (b) for Vela X-1 and GX 1+4. In the panels, the 90% confidence levels are indicated by horizontal dashed lines. In addition, the energies of the iron K-edge for several cases of the ionization state are expressed by vertical dashed lines, calculated from the energy of the iron K_α line at the deepest edge phase and the energy ratio between the absorption edge and the K_α line of iron (Kallman et al. 2004; Yamaguchi et al. 2014). The energy of edge₂ is constrained to 7.23–7.40 keV at the 90% confidence level for GX 1+4, while for Vela X-1, an upper limit of 7.22 keV at the 90% confidence level is determined. These energy ranges correspond to the ionization states of Fe_VI–XI and <Fe_VII, respectively. The determined ionization state of iron within the Alfvén radius for GX 1+4 is equal to or higher than that obtained from the phase-resolved fitting with the assumption of a single absorption edge component, while in the case of Vela X-1, it is the same ionization state as in the phase-resolved analysis. The lower panels of the figure show the sum of the optical depths of the two edge components, τ₁ + τ₂, as a function of the energy of edge₂. As long as the energy of the edge₂ component is in an acceptable range, these values are consistent within the errors with the maximum depths derived from the phase-resolved fitting employing the single absorption edge component, indicated by horizontal dashed lines in the panels.

Zoom In Zoom Out Reset image size

Now, we again discuss the physical properties, number density, and distance from the neutron star of the absorbing matter existing within the Alfvén radius with the ionization state obtained from the above additional analysis. The determined restrictions of the ionization state of iron, Fe_VI–XI, for GX 1+4 roughly correspond to the value of the ionization parameter 44.7 < ξ < 56.2 ($1.65\lt \mathrm{log}\xi \lt 1.75$) in the case of an optically thick plasma (Kallman & McCray 1982, their model 4). The absorbing matter within the Alfvén radius therefore should satisfy the following condition,

$$\begin{eqnarray}&&1.6\times {10}^{35}\,{\mathrm{cm}}^{-1}\lt {{nr}}^{2}\lt 2.0\times {10}^{35}\,{\mathrm{cm}}^{-1},\end{eqnarray} \tag{ 12 }$$

for GX 1+4, so as to satisfy the observed ionization state. Based on this condition, combined with the requirement on the radius (Equation (6)), the acceptable region for the absorbing matter within the Alfvén radius can be represented on the r–n plane. It is indicated by the light blue region in panel (b) of Figure 11.

If the accreting matter along the magnetic field lines at the Alfvén radius has a particle density of n = 10¹⁷ cm⁻³, the ionization state of iron can be Fe_VI–XI, and the thickness of the accretion flow is 10^5–6 cm. This accretion matter can cause the modulating optical depth of the absorption edge with the spin period.

If the modulating optical depth of the absorption edge with the pulse phase originates in the asymmetric structure of the accreting matter corotating with the pulsar spin, a variation of the ionization state of the absorbing matter along the line of sight, namely, a modulation of the energy of the edge, can be expected. Also, a possible modulation of the center energy of the iron K_α line caused by the asymmetric accreting matter is expected due to the variation of the flow velocity along the line of sight. Significantly improved sensitivity and energy resolution from future X-ray missions, such as XRISM, will allow further detailed observations of the absorption edge, as well as of the emission lines.