The Modulating Optical Depth of Photoelectric Absorption Edge with Pulse Phase in Accretion-Powered X-ray Pulsars

TL;DR

This study discovers that the optical depth at the iron K-edge in accretion-powered X-ray pulsars varies with pulse phase, revealing insights into the geometry and ionization state of accreting matter near neutron stars.

Contribution

First detection of pulse phase modulating optical depths at the iron K-edge in X-ray pulsars, linking absorption features to accretion geometry and matter ionization.

Findings

Significant optical depth modulation observed in Vela X-1 and GX 1+4.

Optical depth peaks when X-ray continuum dims.

No significant variation in iron ionization state with pulse phase.

Abstract

We report the first discovery of pulse phase modulating optical depths at the iron K-edge in accretion-powered X-ray pulsars, from Suzaku observations. A significant modulating optical depth of the iron K-edge is detected for Vela X-1 and GX 1+4. Similar trends are seen in GX 301-2 and OAO 1657-415, though with poor statistical significance. The observed iron K-edge exhibits a maximum optical depth, when the X-ray continuum dims, and there is no significant pulse phase variation in the ionization state of iron. The revealed changes in the optical depth with pulse phase can be explained as being due to the accreting matter captured by the magnetic field lines of the pulsar, which co-rotates with the neutron star spin and is responsible for photoelectric absorption. Based on the above interpretation, we propose that the accreting matter within the Alfv\'en radius contains iron with an ionization state of Fe$_{\rm V\hspace{-0.1em}I-X\hspace{-0.1em}I}$ with a particle density of $10^{17}$ cm$^{-3}$ and has a shape flattened along the azimuthal direction, such as an accretion curtain.

Tables

Table 1. Table 1: Log of the Suzaku observations of the four X-ray pulsars.

SourceaaSuzaku studies: GX 301-2 (Suchy et al., 2012), Vela X-1 (Doroshenko et al., 2011; Odaka et al., 2013; Maitra & Paul, 2013), OAO 1657-415 (Pradhan et al., 2014; Jaisawal & Naik, 2014), GX 1+4 (Yoshida et al., 2017)	ObsID	UT start/end time	DurationbbObservational duration time.		ExposureccEffective exposure times with each XIS.				AimddNominal pointing position (XIS or HXD) in the observations (see Mitsuda et al., 2007).	XIS OperationeeClock mode and window option of XIS (see Koyama et al., 2007).
					XIS 0 / XIS 1 / XIS 3					clock,window
		(YYYY-MM-DD hh:mm)	(ks)		(ks)
GX 301-2	403044020	2009-01-05 10:33/2009-01-07 01:00	138\@alignment@align.4	61.8	61\@alignment@align.8	61.8  /	HXD\@alignment@align	normal,1/4  /
Vela X-1	403045010	2008-06-17 04:45/2008-06-18 21:42	147\@alignment@align.4	104.7	104\@alignment@align.7	104.7  /	XIS\@alignment@align	normal,1/4  /
OAO 1657-415	406011010	2011-09-26 09:34/2011-09-28 16:00	195\@alignment@align.9	84.8	84\@alignment@align.8	84.8  /	XIS\@alignment@align	normal,1/4  /
GX 1+4	405077010	2010-10-02 06:43/2010-10-04 12:20	193\@alignment@align.0	97.3	99\@alignment@align.7	88.3  /	HXD\@alignment@align	normal,1/4  /

Equations

$$\displaystyle M(E)=\left\{\begin{array}[]{ll}1&(E<E_{\rm c})\\ \exp\left(-\tau_{\rm edge}\left(\frac{E}{E_{\rm c}}\right)^{-3}\right)&(E\geqslant E_{\rm c}),\end{array}\right.$$

$$R(E)=$$

$$F=\frac{{\chi }_1^2/{\nu }_1}{{\chi }_2^2/{\nu }_2},$$

$$R_{\mathrm{A}}$$

$$r⩽8.1\times 10^8\mathrm{cm}$$

$$r⩽1.4\times 10^9\mathrm{cm}$$

$$n{r}^2>6.9\times 10^{34}{\mathrm{cm}}^{-1},$$

$$n{r}^2>2.0\times 10^{35}{\mathrm{cm}}^{-1}$$

$$\dot{M}$$

$$n(r)$$

$$s(r)=\zeta d_{\mathrm{A}}{R_{\mathrm{A}}}^{1-\gamma }{r}^{\gamma }.$$

$$1.6\times 10^{35}{\mathrm{cm}}^{-1}<n{r}^2<2.0\times 10^{35}{\mathrm{cm}}^{-1}$$

Full text

stars: neutron – X-rays: binaries – pulsars: individual (GX 301-2, Vela X-1, OAO 1657-415, GX 1+4)

1 Introduction

The X-ray emission lines and photoelectric absorption edges, especially of iron, that are often observed in the spectra of accretion-powered X-ray pulsars (e.g., White et al., 1983), can be attributed to the X-ray source and its vicinity and are reflective of the physical properties of the plasma. Therefore, measuring these spectral features can provide information on the physical state of the X-ray source, its vicinity and surrounding matter. The iron emission line at 6.4 keV from accretion-powered X-ray pulsars has been considered to be X-ray fluorescence emitted from neutral or low-ionized iron illuminated by continuous X-rays from the neutron star or its vicinity (e.g., Koyama, 1985; Makishima, 1986). The absorption edge of iron has also been detected in association with the emission line, and the origin of this has been considered to be the same as the fluorescence lines based on the ionization state estimated from the energies of the edge (Makishima, 1986; Ebisawa et al., 1996; Watanabe et al., 2006). It is generally difficult to determine an emission source or an emission region of emission lines, though we can restrict the location of matter responsible for the absorption edge to a region between the X-ray source and the observer. There are other interpretations, however, such as reflective processes. Thus, if the origin of the absorption edge is clarified, the matter distributed around the neutron star can be identified more definitely than a determination based only on the information obtained from the emission lines. A detailed investigation of the absorption edge feature is therefore an interesting tool for studying accretion flow around neutron stars.

The time variation of these spectral features gives various information, such as the size or geometry of the reprocessing region. In accretion-powered X-ray pulsars, the emission line flux sometimes exhibits variability depending on the pulse phase (e.g., Leahy & Matsuoka, 1990; Day et al., 1993; Choi et al., 1996; Vasco et al., 2013). However, sometimes no modulation with the pulse period is found (e.g., Ohashi et al., 1984; Paul et al., 2002; Lei et al., 2009; Suchy et al., 2012). On the other hand, modulation of the optical depth of the absorption edge with the pulse period in X-ray pulsars has never been reported. Pulse phase variations of the absorption column density, which are derived from the low energy cutoff, have been seen in a few accretion powered X-ray pulsars (KS 1947+300;Ballhausen et al. 2016; GX 1+4; Galloway et al. 2001).

Up until the end of its observational operation in 2015, Suzaku observed more than 30 accretion-powered X-ray pulsars, 60 times or more in total. All of the acquired data are available in an archive. From this data we selected pulsars with a spin period longer than 30 s for a phase-resolved analysis. We then selected statistically good datasets, excluding data sets with short exposure time less than 10 ks and faint sources with less than 0.5 count s*-1* on each X-ray imaging spectrometer (XISs; Koyama et al., 2007), giving about 30 available data sets. We extracted the time-averaged spectra from the data. Finally we selected four data sets of four accretion-powered X-ray pulsars, GX 301-2, Vela X-1, OAO 1657-415, and GX 1+4, because their spectra distinctly exhibit the K-shell emission line and absorption K-edge of iron around 6.4 keV and 7.1 keV even without spectral model fitting, as shown in Figure 5.

The Suzaku observations of the four pulsars provide the opportunity to investigate the time variability in both the iron line emissions and the iron absorption edge through phase-resolved spectroscopy thanks to its high sensitivity and excellent energy resolution, as well as their long spin periods. In this paper, we report the results of the detailed spectral analysis of the Suzaku observations. We present the first discovery of modulation in the optical depth of the iron K-edge with the pulsar spin period. These results reveal the physical conditions of the matter surrounding the neutron stars in the systems.

2 Observation and Data Reduction

2.1 Observation

Table 5 summarizes the details of the Suzaku observations of the four selected X-ray pulsars. In all the observations, the XISs were operated in the normal mode incorporating a 1/4 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 or FI CCDs) and 6 keV equivalent electrons for XIS 1 (back-illuminated or 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 $\mu$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 this data because of its poor statistics for phase-resolved analysis with a spin period of 700 s.

2.2 Data Reduction

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 broad-band 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 (version 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 3′radius and an annulus with an inner radius of 4′and an outer radius 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 0.′3 and 0.′9, respectively, to eliminate the pileup effect. This procedure reduces the pileup fraction to below 1%. The redistribution matrix files 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.

3 Analysis and Results

3.1 Light Curves and Pulse Period

For the timing analysis, we applied a barycentric correction to the arrival times of individual photons using the aebarycen task of FTOOLS (Terada et al., 2008). As for HXD events, dead-time corrections were applied using an 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., 1/4 window option). Data from all 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 5 shows 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 $\sim 1.4\times 10^{5}$ s and a high intensity period followed for $\sim 6\times 10^{4}$ 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\times 10^{5}$ 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 $\sim 5\times 10^{4}$ seconds of the observation, when the source exhibited high intensity. The selected time interval is shown in Figure 5 by the shaded area.

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 LABEL:tab:TimAveResult. Errors estimated according to Leahy (1987) are also given in the table. All the sources are consistent with the previous reported values derived from the same Suzaku data (Suchy et al., 2012; Maitra & Paul, 2013; Pradhan et al., 2014; Yoshida et al., 2017), within the errors.

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 keV and 7.1–12.0 keV ranges are shown in panels (a) and (b) of Figure 3–6, where the data from the three XISs are summed. According to the criteria described in § 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.

3.2 Time-averaged Spectrum

Figure 5 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 § 3.1. The data of Vela X-1 obtained by the three XISs was 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_{\rm H}=10^{23}$ H atoms cm*-2*), as seen in Figure 5, 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 is not well calibrated (Suchy et al., 2012). Data of XIS in the 2.2–2.4 keV energy range was 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 fitted by the model employed by earlier studies, to describe the same Suzaku broad-band spectra (Suchy et al., 2012; Pradhan et al., 2014; Maitra & Paul, 2013; Yoshida et al., 2017). The spectral continuous components used are summarized in Table LABEL:tab:TimAveResult. 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 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_{\rm Fe}$, to vary, and values of $Z_{\rm Fe}$ $\sim 12$ in Vela X-1 and $\sim 1.4$ in GX 1+4 were typically obtained. The obtained large abundance of Vela X-1 will be discussed in 3.3.1.

The results of the time-averaged spectral fitting are summarized in Table LABEL:tab:TimAveResult. 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 LABEL:tab:TimAveResult. The distances to the sources from the observer are also given in the table. These values are consistent with those reported by Suchy et al. (2012), Doroshenko et al. (2011), Jaisawal & Naik (2014), and Yoshida et al. (2017) for the same Suzaku data.

3.3 Phase-resolved Analysis

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 the figure) were combined to improve the photon statistics.

We fitted the broad-band phase-resolved spectra with the same models as those used in the time-averaged spectral fitting listed in Table LABEL:tab:TimAveResult. 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_{\rm Fe}$ was fixed at the value obtained from the time-averaged fitting. We found that the adopted model fitted all 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, 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 parameter. 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 ${\rm K_{\alpha}}$ 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$\pm$2 of Vela X-1 obtained from the time-averaged spectral fit, shown in Table LABEL:tab:TimAveResult, 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 broad-band 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 and 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:

[TABLE]

where $E_{\rm c}$ is the threshold energy and $\tau_{\rm edge}$ is the absorption optical depth at the threshold energy. The applied model consists of ${\tt powerlaw}\times{\tt edge}+{\tt gaussian}+{\tt gaussian}+{\tt gaussian}$, in which three gaussian models express the iron ${\rm K_{\alpha}}$, iron ${\rm K_{\beta}}$, and nickel ${\rm K_{\alpha}}$ emission lines. Since the iron ${\rm K_{\beta}}$ 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 ${\rm K_{\alpha}}$ emission line is around 6.4 keV, the ionization state of the iron should be low (FeXVII at most). In the above ionization state, the difference of the ratio between the energies of the iron ${\rm K_{\alpha}}$ and iron ${\rm K_{\beta}}$ lines is at most 1%, but that of the energy of the iron edge to that of iron ${\rm K_{\alpha}}$ line is at most 7% (Kallman et al., 2004; Yamaguchi et al., 2014). Therefore, we fixed the ratio of the energy of the iron ${\rm K_{\beta}}$ line to that of the iron ${\rm K_{\alpha}}$ 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 ${\rm K_{\alpha}}$ 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 ${\rm K_{\alpha}}$ 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 broad-band 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 the CCDs.

Parameters obtained from the spectral fitting are summarized in Table LABEL:tab:phasedivfit and plotted as a function of pulse phase in Figures 3–6. We find that the optical depths of the iron K-edge, $\tau_{\rm edge}$, exhibit modulations with pulse phase, as seen in panel (e) of each figure. For all the four pulsars, $\tau_{\rm edge}$ increases at the phase when the X-ray continuum flux becomes dim. The energies of the iron ${\rm K_{\alpha}}$ 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 ${\rm K_{\alpha}}$ 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 $\tau_{\rm edge}$. The expected energy ratios for FeI–FeVII 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${}_{\rm II\mbox{\textendash}VII}$. We confirmed that the assumed ratio of the energy of the iron ${\rm K_{\beta}}$ line to that of the iron ${\rm K_{\alpha}}$ line of 1.103 in the fitting is consistent with this ionization state. Table LABEL:tab:phasedivfit lists the obtained intensities of the iron ${\rm K_{\alpha}}$ and ${\rm K_{\beta}}$. The ratio 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 of interval 3 of OAO 1657-415. In these cases, the assumed energy ratio between the iron ${\rm K_{\alpha}}$ and ${\rm K_{\beta}}$ lines might be smaller than the actual value, and 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 ${\rm K_{\beta}}$ line.

Using the obtained optical depth, an equivalent hydrogen column density for an assumed solar abundance of $N_{\rm H}=\frac{\tau_{\rm edge}}{Z_{\rm Fe}\sigma_{\rm Fe}}$ is derived, where $Z_{\rm Fe}$ and $\sigma_{\rm Fe}$ are the iron abundance in the solar system ($=3.2\times 10^{-5}$; Asplund et al., 2009) and the photoelectric absorption cross section of FeII at the energy of the K-edge ($=3.3\times 10^{-20}$ ${\rm cm}^{2}$; Verner & Yakovlev, 1995). The estimated equivalent hydrogen column densities for each phase bin are given in Table LABEL:tab:phasedivfit.

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 $\chi^{2}$ tests to evaluate the statistical significance, with a null hypothesis that there is no variation of the optical depth over the pulse phase. The $\chi^{2}$ values obtained from the constant fittings to the optical depth in the four pulsars are summarized in Table LABEL:tab:sumstattest 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 keV and 7.5 keV, resulting from variations between the two phases of the equivalent widths of the iron ${\rm K_{\alpha}}$ line and nickel ${\rm K_{\alpha}}$ line.

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

[TABLE]

where $A$ is the slope, $B$ is the intercept, $M(E)$ is the edge component given by Equation 3, and $G_{\rm Fe}(E)$ and $G_{\rm Ni}(E)$ are Gaussian functions expressing the emission line features at 6.4 keV, and 7.5 keV. In this case, the $\tau_{\rm edge}$ in $M(E)$ roughly corresponds to the difference in the optical depth at the iron K-edge between the two spectra, denoted by $\Delta$. 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_{\rm Ni}(E)$ to that of $G_{\rm Fe}(E)$ is fixed at 1.168 for the neutral case, and the width of $G_{\rm Ni}(E)$ is assumed to be equal to that of $G_{\rm 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 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 $\Delta$ derived from the fitting are listed in Table LABEL:tab:sumstattest. It is confirmed that the values of $\Delta$ are consistent with the differences between the optical depths of the two phase-resolved spectra, within the errors, for all sources (see Table LABEL:tab:phasedivfit). The resultant $\chi^{2}$ values are also summarized in Table LABEL:tab:sumstattest.

To evaluate the statistical significance of the $\chi^{2}$ 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

[TABLE]

where $\chi_{1}^{2}$ and $\chi_{2}^{2}$ are chi-squared values and $\nu_{1}$ and $\nu_{2}$ 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 LABEL:tab:sumstattest. The corresponding chance probabilities in the random case of an improvement of $\chi^{2}$, given in Table LABEL:tab:sumstattest, are not small, and we can only marginally conclude the differences of the optical depths by the F-test.

However, the residuals shown in the 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_{\rm 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 reproduce only 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_{\rm 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_{\rm 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 are summarized in Table LABEL:tab:sumstattest, 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, and 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 significances from the run-test.

4 Discussion

Focusing on the absorption edges and the emission lines of iron in X-ray pulsars spectra, we conducted pulse phase-resolved spectroscopy with data derived from Suzaku observations. We found a significant modulating optical depth of the iron K-edge with the pulse phase from Vela X-1 and GX 1+4. Similar trends were seen for GX 301-2 and OAO 1657-415, though with poor statistical significance. Although the relation between this modulation and the absorption column density derived from the low energy cutoff is interesting, we found that the obtained absorption column densities strongly correlate with the power-law index and we consider that the obtained values are not physically meaningful. Therefore, this relation is not discussed in this work. In the following, we discuss the origin of the absorbing matter based on the parameters obtained for Vela X-1 and GX 1+4. Hereafter we refer to the two phases at which the observed the iron K-edge exhibits the maximum and minimum optical depth as the deepest and shallowest edge phases, respectively.

We assume that the same matter is responsible for the absorption and fluorescent emission, which is illuminated by X-rays from the vicinity of the neutron star. Then, we can estimate the ionization state of iron from the ratios of the absorption edge energy to the emission line energy. We determined the ionization state to be Fe${}_{\rm II\mbox{\textendash}VII}$ over the pulse phase. This ionization state is common between the two pulsars. From these ionization states, the ionization parameter $\xi\equiv\frac{L_{\rm X}}{nr^{2}}$ (Tarter et al., 1969, where $L_{\rm X}$ is the luminosity of the X-ray source, $n$ is the number density of the gas, and $r$ is the distance between the X-ray source and the gas) can be estimated. The estimated ionization state of iron corresponds to a value of $\xi$ of less than 44.7 ($\log\xi<1.65$) in the case of an optically thick plasma (Kallman & McCray, 1982, their model 4, in which the incident X-rays are assumed to be 10-keV bremsstrahlung with $L_{\rm X}=10^{37}$ erg s*-1*, roughly simulating the GX 1+4 case). In the following discussion, we therefore require the ionization parameter of the absorbing matter, so as to satisfy the observed ionization state. Unless otherwise specified, the ionization parameter is restricted to $\xi<44.7$.

4.1 Possible Origin of Modulating Optical Depth of Absorption Edge

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

The variation in the physical states of the absorbing matter, such as ionization state and density, according to the spin phase, is a possible explanation for the modulating optical depth of the absorption edge with the pulse phase. Since X-rays from the pulsar are not isotropic, the matter located between the observer and the X-ray source is strongly irradiated by the X-rays, when the observed X-ray flux of the pulsar is intense. The irradiated matter might hence change its ionization state. This means that the variation in the optical depth at the absorption edge may result from a change in the number of iron ions in a certain ionization state, which are responsible for the absorption edge of interest. In this interpretation, the deepest edge phase should correspond to the phase for which the observed X-ray flux from the pulsar decreases, while the shallowest edge phase should occur along with an increase of observed X-ray flux from the pulsar. In addition, an enhancement of the illuminating X-ray intensity to the gas in the line of sight increases the ionization state of the gas, and the resulting energy of the absorption edge should become high. On the other hand, the gas recombines as the illuminating X-ray intensity declines, entering a low ionized state, and we can expect the energy of the absorption edge to become low. As a result, this effect can lead to a modulation of the energy with the X-ray flux.

However, no modulating energy of the edge in association with increasing or decreasing continuous X-ray flux was observed with Suzaku. The ratio of the iron K-edge energy to that of the iron ${\rm K_{\alpha}}$ line did not change significantly throughout the pulse phase for Vela X-1 (see Figure 4). Furthermore, the energy ratio of GX 1+4 becomes large at the deepest edge phase, where the continuous flux becomes minimum (see Figure 6). This result is the opposite trend than that predicted for the above effect. These observational results demonstrate that a variation in the physical state of the absorbing matter in the line of sight is not a plausible explanation for the modulating optical depth of the absorption edge with the pulse phase.

4.1.2 Modulating Reflected X-rays

If there is neutral matter surrounding the pulsar, a fraction of the X-rays from the neutron star will be reflected by the neutral matter, e.g. a part of the accretion disk, and a composite X-ray spectrum will be observed (e.g., Lightman & White, 1988). In fact, X-ray reflection on accretion disks is observed in neutron star binaries (Rea et al., 2005; Cackett et al., 2010). The reprocessed X-rays intrinsically exhibit a distinct iron edge feature if the reflecting matter has a composition comparable to the solar composition. In addition, a neutral iron emission line coming from the reprocessing matter can be observed.

We now consider X-ray reflection from neutral matter exposed to X-rays from an X-ray pulsar, whose emission profile is not spherically symmetric. In this case, the appearance of the reflected X-ray component can change according to the rotation of the neutron star. As a consequence, this modulation of the reflected X-rays can cause variations with the pulse period in the optical depth of the iron edge as well as the flux of the iron emission line of the observed X-ray spectrum. It is expected that when an increasing X-ray reflection component is observed, the optical depth of the edge and the emission line flux will be simultaneously enhanced, in a simple case.

For GX 1+4, although modulations of both the optical depth of the edge and the emission line flux (Yoshida et al., 2017) are observed in the Suzaku observation, their peak-phases are not synchronized with each other. In the case of Vela X-1, a variation in the optical depth of the edge is detected with statistical significance, while the measured line flux did not change significantly throughout the pulse phase.

However, since there might be other causes changing the emission line flux, we cannot judge from these observational results whether the modulating reflected X-rays can be attributed to the modulating optical depth of the edge with the pulse phase.

4.1.3 Movement of the Emission Region

A change in the amount of matter along the line of sight by a movement of the emission region or of the absorbing matter itself, due to the neutron star spin, can possibly explain the modulating optical depth of the absorption edge with the pulse phase. If the matter is moving with the spin of the neutron star, the matter should exist at or inside the Alfvén radius (Lamb et al., 1973; Elsner & Lamb, 1977). This case will be discussed in the next section.

In this section, we suppose that the matter does not move according to the neutron star spin, and hence the absorbing matter exists outside the Alfvén radius. A modulating optical depth with the pulse phase may appear if the emission region moves according to the neutron star spin. The emission region is thought to exist around the magnetic poles of a neutron star and its size is considered to be less than or comparable to the neutron star radius, namely $10^{6}$ cm. On the other hand, the distance of the absorbing matter from the emission region is expected to be more than the Alfvén radius, namely $10^{8}$ cm. By considering this geometry, even if the position of the emission region changes, no change of the absorption column density with the pulse phase can be expected, unless there is a very stable and tiny structure in the absorbing matter, which exists outside the Alfvén radius.

4.1.4 Change of Absorbing Matter Location

Alternatively, if the absorbing matter co-rotates with the neutron star spin, variation of the optical depth with the pulse phase can occur. The co-rotation with the neutron star spin implies that the absorbing matter should be confined by the magnetic field of the neutron star, and hence it should be located at or within the Alfvén radius $R_{\rm 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

[TABLE]

where $M_{\ast}$, $R_{6}$, $B_{12}$, and $L_{37}$ are the mass, radius (denoted as $R_{\ast}$) in units of $10^{6}$ cm, surface magnetic field strength (denoted as $B_{\rm surf}$) in units of $10^{12}$ G of the neutron star, and luminosity (denoted as $L_{\rm X}$) in units of $10^{37}$ erg s*-1*, respectively (Frank et al., 2002). Substituting a neutron star radius of $R_{\ast}=10$ km and a mass of $M_{\ast}=1.4M_{\odot}$ in Equation 6, 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\times 10^{12}$ 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^{13}$ G, since no cyclotron resonance scattering features have been reported in its spectra (Yoshida et al., 2017).

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

[TABLE]

for Vela X-1 ($L_{\rm X}=3.1\times 10^{36}\ {\rm erg\,s^{-1}}$ and $B_{\rm surf}=2.2\times 10^{12}\,{\rm G}$) and

[TABLE]

for GX 1+4 ($L_{\rm X}=8.9\times 10^{36}\ {\rm erg\,s^{-1}}$ and $B_{\rm surf}=1.0\times 10^{13}\,{\rm 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

[TABLE]

and

[TABLE]

for Vela X-1 and GX 1+4, respectively, where $n$ is the number density of particles of the matter. The acceptable regions on the $r$–$n$ planes satisfying the requirement on the absorbing matter given by the above equations, 7, 8, 9, and 10, 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 $\Delta N_{\rm H}=10^{23}$ cm*-2*, the geometric thickness of the absorbing matter along the line of sight, $\delta$, 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^{18}$ ${\rm cm}^{-3}$, the absorbing matter contains almost neutral irons and forms a structure whose geometric thickness along the line of sight is $10^{5}$ cm. Given these facts, the observed variation in the amount of absorption in the line of sight due to the absorbing matter co-rotating with the neutron star spin is a plausible explanation for the modulating optical depth of the iron absorption edge with the pulse phase.

4.2 Accretion Flow as Absorbing Matter

The condition described in § 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

[TABLE]

where $s(r)$ is the cross sectional area of the accretion flow, $m_{\rm p}$ is the proton mass, and $\mu_{\rm 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 free fall velocity given by $v(r)=\sqrt{2GM_{\ast}/r}$. By assuming that all the kinetic energy liberated is converted to radiation energy at the neutron star surface, namely $L_{\rm X}=GM_{\ast}\dot{M}/R_{\ast}$, the density in the accretion flow is given by

[TABLE]

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^{\gamma}$. A schematic picture of the assumed accretion flow is shown in Figure 12. We assume that the width is $w(r)=\zeta r$, where $\zeta$ 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_{\rm A}$ is the thickness at $R_{\rm A}$. Then, we obtain

[TABLE]

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 $\zeta=0.4\pi$ (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 $\gamma$ 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 $\gamma=3$. Since this assumption of $\gamma=3$ can be applied only near the neutron star and cannot be extrapolated to the Alfvén radius, we consider that the actual $\gamma$ value should be between 2 and 3. Thus, we calculated both cases of $\gamma=2$ and 3. Substituting Equation 13 into Equation 12 with $\gamma=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 $\zeta=0.4\pi$, the calculated number densities in the accretion flow are plotted by red dashed lines ($\gamma=2$) and blue dotted lines ($\gamma=3$) in Figure 11 for some typical thicknesses for the accretion flow at the Alfvén radius, $d_{\rm A}=10^{4}$, $10^{5}$, and $10^{6}$ 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_{\rm A}<10^{5}$ cm and its cross sectional area is proportional to $r^{3}$. On the other hand, for $\gamma=2$, only the outer region is acceptable for accreting matter of thickness $d_{\rm A}\sim 10^{4\mbox{\textendash}5}$ cm, for both sources. It should be noted that the thickness of the modeled accretion flow at the Alfvén radius, $d_{\rm A}$, is roughly consistent, within a factor of 3–8, with $\delta$, indicated on the vertical axes on the right of the panels in Figure 11.

As mentioned in § 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; Giles et al., 2000; Galloway et al., 2000, 2001). 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 13, its thickness is sufficiently smaller than its width if the thickness is the estimated value of $10^{4{\mbox{\textendash}}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, latter shows the schematic figure of accretion geometry) may be more plausible than a whole cone geometry.

4.3 Distribution of Iron Surrounding X-ray Pulsar

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 co-rotates 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_{\rm H}=10^{23}$ cm*-2*, was observed (see Table LABEL:tab:phasedivfit). 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\times 10^{11}$ 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 $\sim$30 s (Yoshida et al., 2017). A possible explanation for this delay was proposed by Yoshida et al. (2017), where the iron line is a fluorescence of the matter exposed to the X-ray continuum from the neutron star, at a distance of $\sim 10^{12}$ 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.

4.4 Two Absorption Edges Hypothesis

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 § 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 ${\rm K_{\alpha}}$ 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 § 3.3. For the spectral fitting in the restricted energy range, a model consisting of ${\tt powerlaw}\times{\tt edge_{1}}\times{\tt edge_{2}}+{\tt gaussian}+{\tt gaussian}+{\tt gaussian}$ is applied to express the absorption by iron in the two regions. In this model, ${\tt edge_{1}}$ and ${\tt edge_{2}}$ 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 ${\rm K_{\beta}}$ line to that of the iron ${\rm K_{\alpha}}$ line was fixed to 1.103 for the neutral case (Yamaguchi et al., 2014). The center energy of the nickel ${\rm K_{\alpha}}$ line was also fixed to a value derived from the phase-resolved fitting at the deepest edge phase. The energy and the optical depth of the ${\tt edge_{1}}$ component were fixed to the values obtained from the fitting of the shallowest edge phase. The energy of the ${\tt edge_{2}}$ 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 ${\tt edge_{2}}$ component was allowed to be free. The resultant $\chi^{2}$ values as a function of the energy of the ${\tt edge_{2}}$ are plotted in the upper panels of Figure 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 ${\rm K_{\alpha}}$ line at the deepest edge phase and the energy ratio between the absorption edge and the ${\rm K_{\alpha}}$ line of iron (Kallman et al., 2004; Yamaguchi et al., 2014). The energy of ${\tt edge_{2}}$ are 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${}_{\rm VI{\mbox{\textendash}}XI}$, and $<$FeVII, 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, $\tau_{1}+\tau_{2}$, as a function of the energy of ${\tt edge_{2}}$. As long as the energy of the ${\tt edge_{2}}$ 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.

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${}_{\rm VI{\mbox{\textendash}}XI}$, for GX 1+4 roughly correspond to the value of the ionization parameter $44.7<\xi<56.2$ ($1.65<\log\xi<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 conditions:

[TABLE]

for GX 1+4, so as to satisfy the observed ionization state. Based on this condition, combined with the requirement on the radius (Equation 8), 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^{17}$ cm*-3*, the ionization state of iron can be Fe${}_{\rm VI{\mbox{\textendash}}XI}$ and the thickness of the accretion flow is $10^{5{\mbox{\textendash}}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 co-rotating 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 ${\rm K_{\alpha}}$ 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.

5 Conclusion

We performed a study of X-ray pulsars focusing on the absorption edge and the emission line of iron, and investigated the properties of the iron surrounding the pulsar and the accretion flow along the magnetic field.

We discovered a significantly modulating optical depth of the iron K-edge with the neutron star spin for Vela X-1 and GX 1+4. Similar trends were found for GX 301-2 and OAO 1657-415, though with poor statistical significance. A possible interpretation of the observed changes of the optical depth with pulse phase is that the accreting matter captured by the magnetic field lines of the pulsar, which co-rotates with the neutron star spin, causes photoelectric absorption. Based on the observations of the iron absorption edge, we speculate that iron surrounding an X-ray pulsar is distributed in two regions divided by the Alfvén radius, where the iron is in different physical conditions of ionization state and geometry. In particular, we propose that the accretion flow along the magnetic field line of a neutron star within the Alfvén radius contains iron with an ionization state of Fe${}_{\rm VI{\mbox{\textendash}}XI}$, with a particle density of $n\sim 10^{17}$ cm*-3*.

This research was carried out using data obtained from the Data Archive and Transmission System (DARTS), provided by the Center for Science-satellite Operation and Data Archive (C-SODA) at ISAS/JAXA. This work was partially supported by a MEXT-Supported Program for the Strategic Research Foundation at Private Universities, 2014–2018 (S1411024).