A CLUMPY STELLAR WIND AND LUMINOSITY-DEPENDENT CYCLOTRON LINE REVEALED BY THE FIRST SUZAKU OBSERVATION OF THE HIGH-MASS X-RAY BINARY 4U 1538−522

To mitigate pile-up in the detectors, the three XIS instruments operated in 1/4 window mode. Each XIS unit features ⁵⁵Fe calibration sources illuminating two corners of each CCD, but this windowing mode precluded their use. An additional attitude correction was performed on the event files using the Suzaku tool aeattcor2 prior to spectral extraction. A pileup estimation using the pileest tool showed a maximum pileup of ∼7% around the center of the image, and so the source regions for each XIS unit were defined as annuli, excising the central region down to ∼3% pileup. Background regions were defined as strips along the edges of each image; in the case of XIS0, the background regions additionally avoided a strip of bad pixels located along one edge of the image (due to what was likely a micrometeorite impact in 2009—see Tsujimoto et al. 2010). In all detectors, background regions were defined with the same total area as the source region. Source and background spectra and light curves for each XIS unit were extracted using xselect, totaling 46 ks of exposure in each detector, and response matrices and ancillary response files were generated using xisrmfgen and xissimarfgen. Source spectra were regrouped using grppha; the grouping used is the same as that used by Nowak et al. (2011), which attempts to optimally account for the energy resolution of the XIS CCDs. Data below 1 keV and between 1.6 and 2.3 keV showed significant discrepancies between XIS1 and XIS0+XIS3, and were ignored; we additionally ignored XIS data above 10 keV.

Data from the Suzaku HXD/PIN were extracted using xselect. The "tuned" non-X-ray background (NXB) event files provided by the Suzaku team were used, while the cosmic X-ray background (CXB) was simulated according to the model of Boldt (1987), using the fakeit procedure in XSPEC v.12.8.0. The NXB and CXB spectra were added using mathpha to produce a single PIN background spectrum. The source spectrum was deadtime-corrected using the hxddtcor tool; the total PIN exposure after deadtime correction is 34 ks. PIN spectra were regrouped to 100 counts/bin, and data were ignored below 15 keV in all data sets and above ∼40–50 keV, depending on the signal-to-noise in the particular data set being analyzed.

We do not examine the full observation in a single spectrum due to the high variability seen with the onset of the dip. Rather, we investigate here the time-averaged spectrum of the first 32 ks of the observation, where the source's spectral variability is relatively low (see Sections 4.2 and 4.3). The three XIS and single PIN spectra cover an energy range of 1 to 40 keV, with gaps between 1.6 and 2.3 keV (due to uncertainties in the instrumental response) and between 10 and 15 keV (where both the XIS and the PIN lack good coverage). We fit our broadband spectra with three different continuum models: a single power law multiplied by either the highecut (White et al. 1983) or fdcut (Tanaka 1986) exponential cutoff, and the Negative–Positive EXponential model (npex; Mihara 1995), with two power laws:

$$\begin{equation} \mathtt {pow*highecut}(E) = \left\lbrace \begin{array}{@{}ll}A E^{-\Gamma } & E < E_{\mathrm{cut}} \\ A E^{-\Gamma }\exp \left(\frac{E_{\mathrm{cut}} - E}{E_{\mathrm{fold}}} \right) & E \ge E_{\mathrm{cut}} \\ \end{array}\right.\quad \end{equation} \tag{ 1 }$$

$$\begin{equation} \mathtt {pow*fdcut}(E) = A E^{-\Gamma } \left[1 + \exp \left(\frac{E - E_{\mathrm{cut}}}{E_{\mathrm{fold}}}\right)\right]^{-1} \end{equation} \tag{ 2 }$$

$$\begin{equation} \mathtt {npex}(E) = A (E^{-\alpha } + BE^{+\beta }) e^{E/kT}. \end{equation} \tag{ 3 }$$

The remaining widely used power-law/exponential cutoff continuum model, cutoffpl, did not provide good fits, with significant additional structures at higher energies and a reduced χ² of ∼2, and so it was not used. In Equations (1) and (2), the parameter Γ is the power-law index. We should note here that despite using the same parameter names, the E_cut and E_fold parameters of the highecut and fdcut models are mathematically different and should not, in principle, be compared directly (although they do play similar roles). npex (Equation (3)) has two power laws, with the α index being negative-definite and the β index being positive-definite. The B parameter is the relative normalization of the positive index power law. The considerable added leverage of the second power law tends to make fitting with npex a difficult proposition; we avoid this somewhat by freezing β to 2.0; this results in the model simulating the Wien hump at E ∼ kT that should exist in the Comptonized spectrum of the neutron star (Makishima et al. 1999). npex features a kT parameter, which functions in a similar capacity to E_fold in the other two models.

Each continuum model is modified by absorption, using the latest version ⁸ of the tbnew absorption model (Wilms et al. 2010) using the abundances of Wilms et al. (2000) and cross sections from Verner et al. (1996). In addition to the broad effects of absorption, the residuals have significant additional structures between 6.4 and ∼7.5 keV and at ∼23 keV. The strong emission line at ∼6.4 keV is the Kα emission of neutral iron, while additional emission and absorption-like features in the residuals at ∼7 keV are identified as the iron Kβ emission line (7.06 keV) and the iron K-shell ionization edge (7.112 keV in neutral iron). The emission lines are modeled with additive Gaussians while the edge feature is modeled by allowing the abundance of iron in the absorber to vary. We fix the widths of the emission lines to 0.01 keV, as the line widths could not be constrained. We could not resolve the 6.4 keV feature into multiple emission lines as was done by Rodes-Roca et al. (2010). At ∼23 keV there is a broad absorption-like feature, which we identify as 4U 1538−522's CRSF. This feature is included in the model using a local XSPEC model, gauabs, a multiplicative absorption feature with a Gaussian optical depth profile:

$$\begin{equation} \mathtt {gauabs}(\tau) = e^{-\tau \left(E\right)} \end{equation} \tag{ 4 }$$

$$\begin{equation} \tau (E) = \tau _{0}\exp \left(-\frac{\left(E-E_{0}\right)^{2}}{2\sigma ^{2}}\right). \end{equation} \tag{ 5 }$$

Here, E₀ is the centroid energy of the feature, τ₀ is the maximum optical depth, and σ is the width of the optical depth profile. The gauabs component modifies only the power law and exponential cutoff continuum, as it is produced via the same physical processes that produce the rest of the continuum deep in the accretion column. The two Gaussian emission lines, Fe Kα and Kβ, are added to this continuum, and the result is then modified by the tbnew absorption model. Unabsorbed fluxes are determined using the cflux model in XSPEC. The full model additionally has a multiplicative constant in order to account for calibration differences between instruments. The normalization constant for XIS0 was frozen at 1, while the constants for XIS1 and XIS3 were allowed to vary. The PIN normalization was frozen at the nominal value of 1.16 per the Suzaku ABC guide, as the lack of overlap between the PIN and XIS made this parameter poorly constrained (and influential over values of the high-energy cutoff parameters). The fitted parameters with 90% error bars can be found in Table 1, and the best-fit spectrum with the fdcut continuum is plotted in Figure 3.

Zoom In Zoom Out Reset image size

Table 1. Spectral Fits for Time-averaged Spectrum

Parameter	Units	highecut	fdcut	npex
Fluxa	× 10−10 erg cm−2 s−1	4.2 ± 0.2	4.2 ± 0.2	4.2 ± 0.2
NH	× 1022 cm−2	2.1 ± 0.2	2.1 ± 0.2	1.7 ± 0.3
Γ		1.17 ± 0.01	$1.16^{+ 0.01}_{- 0.02}$	...
Ab	× 10−2 photons cm−2 s−1	4.75 ± 0.08	4.70 ± 0.08	4.76 ± 0.18
Ecut	keV	$21^{+2}_{-3}$	27 ± 1	...
Efold	keV	$8^{+2}_{-1}$	5.0 ± 0.7	...
α		...	...	0.67 ± 0.03
B	× 10−3	...	...	3.2 ± 0.2
kBT	keV	...	...	$4.83^{+ 0.09}_{- 0.08}$
EKα	keV	6.426 ± 0.008	6.426 ± 0.008	6.425 ± 0.008
IKαc	× 10−4 photons cm−2 s−1	3.7 ± 0.3	3.6 ± 0.3	3.7 ± 0.3
EKβ	keV	$7.05^{+0.08}_{-0.07}$	7.05 ± 0.08	$7.06^{+0.13}_{-0.07}$
IKβd	× 10−5 photons cm−2 s−1	4 ± 2	3 ± 2	4 ± 2
Fe abundancee		$1.4^{+ 0.9}_{- 0.7}$	$1.2^{+ 0.9}_{- 0.7}$	>1.06
Ecyc	keV	$22.2^{+ 0.8}_{- 0.7}$	23.0 ± 0.4	22.4 ± 0.3
σcyc	keV	3.0 ± 0.3	3.2 ± 0.4	2.9 ± 0.4
τcyc		$0.9^{+ 0.1}_{- 0.2}$	$0.72^{+ 0.09}_{- 0.08}$	0.58 ± 0.05
$\chi ^{2}_{\mathrm{red}}$ (dof)		1.17(444)	1.21(444)	1.23(444)

Notes. ^a2–10 keV flux. ^bPower-law normalization at 1 keV. ^cIron Kα intensity. ^dIron Kβ intensity. ^eRelative to ISM abundances.

Download table as:  ASCIITypeset image

Gaps in the light curve, caused by the passage of the satellite through the South Atlantic Anomaly (SAA), provide convenient points for a time-resolved spectral analysis. We thus divide the observation into nine segments, defined by the SAA gaps. The exception to this is for the fifth and sixth time bins, where the dividing line is placed at the start of the flare rather than during the SAA gap. This ensures that the fifth bin only encompasses the dip, and the sixth, the flare. The time intervals for these spectra are indicated in Figure 1. The HXD/PIN and three XIS spectra for each time bin were initially fit with the same set of models used for the time-averaged spectrum. The iron Kβ was generally not significantly detected in these spectra, but we did include it in the model (with its energy fixed at 7.056 keV) in order to obtain upper limits on its intensity. The iron abundance was additionally fixed to the ISM abundance, as it was unconstrained at the high end.

All three continuum models fit the pre-dip and spike spectra well, with $\chi ^{2}_{\rm red} \sim 1.0\hbox{--}1.1$. The dip and the final two spectra (spectra 4, 8, and 9 in Table 2), however, show noticeable curvature in their XIS spectra and have much poorer fits, with the dip spectrum having a reduced χ² of 1.6 for 427 degrees of freedom and spectra 8 and 9 having reduced χ² of ∼1.3 for 415 and 421 degrees of freedom, respectively. We thus modify the model for these spectra by applying a partial-covering absorption component of the form

$$\begin{equation} \mathtt {partcov}(E) = (1-f_{{\rm pcf}})*\mathtt {tbnew}_{1}(E) + f_{{\rm pcf}}*\mathtt {tbnew}_{2}(E). \end{equation} \tag{ 6 }$$

Table 2. Spectral Fits for Time-resolved Spectroscopy

Parameter	Units	1	2	3	4	5	6	7	8	9
highecut
Fluxa	× 10−10 erg cm−2 s−1	4.29 ± 0.03	5.48 ± 0.05	4.59 ± 0.05	4.20 ± 0.05	4.4 ± 0.2	8.7 ± 0.1	4.65 ± 0.05	$2.79^{+ 0.10}_{- 0.09}$	3.6 ± 0.1
NH	× 1022 cm−2	2.22 ± 0.03	2.14 ± 0.04	2.27 ± 0.05	1.98 ± 0.05	3.3 ± 0.2	2.82 ± 0.08	2.16 ± 0.05	2.6 ± 0.2	3.1 ± 0.2
NH, pcf	× 1022 cm−2	...	...	...	...	17 ± 3	...	...	$8^{+4}_{-2}$	$9^{+4}_{-3}$
fpcf		...	...	...	...	0.61 ± 0.03	...	...	0.40 ± 0.09	$0.41^{+ 0.09}_{- 0.08}$
Γ		1.16 ± 0.01	1.18 ± 0.02	1.21 ± 0.02	1.20 ± 0.02	1.13 ± 0.06	0.94 ± 0.03	1.11 ± 0.02	$1.38^{+ 0.06}_{- 0.07}$	1.38 ± 0.06
Ab	× 10−2 photons cm−2 s−1	4.4 ± 0.1	5.8 ± 0.2	5.1 ± 0.2	4.6 ± 0.2	$4.2^{+ 0.6}_{- 0.5}$	4.4 ± 0.1	4.3 ± 0.2	$4.1^{+ 0.6}_{- 0.5}$	$5.3^{+ 0.7}_{- 0.6}$
Ecut	keV	$24^{+2}_{-4}$	$23^{+3}_{-4}$	$15^{+3}_{-2}$	22 ± 5	$21^{+4}_{-5}$	14.1 ± 0.8	$18^{+8}_{-5}$	$17^{+6}_{-2}$	$17^{+8}_{-2}$
Efold	keV	$6^{+2}_{-1}$	$7^{+3}_{-2}$	$13^{+2}_{-3}$	$9^{+4}_{-3}$	$8^{+4}_{-2}$	12 ± 1	$11^{+4}_{-5}$	$12^{+3}_{-5}$	$12^{+2}_{-6}$
EKα	keV	6.43 ± 0.01	6.43 ± 0.02	6.44 ± 0.02	6.40 ± 0.02	6.41 ± 0.01	6.39 ± 0.02	6.40 ± 0.02	6.42 ± 0.04	6.43 ± 0.03
IKαc	× 10−4 photons cm−2 s−1	$3.6^{+ 0.4}_{- 0.3}$	4.4 ± 0.7	3.8 ± 0.6	3.3 ± 0.6	4.9 ± 0.5	8 ± 1	3.2 ± 0.6	1.2 ± 0.5	1.9 ± 0.5
IKβd	× 10−5 photons cm−2 s−1	4 ± 3	<8.99	<9.39	<7.86	9 ± 5	<15.79	<12.43	<4.12	8 ± 5
Ecyc	keV	23 ± 1	23 ± 1	$21^{+3}_{-1}$	$22^{+2}_{-1}$	$21.7^{+1.0}_{- 0.9}$	23 ± 1	$21.3^{+1.0}_{- 0.7}$	$21.4^{+ 0.7}_{- 0.6}$	$21.8^{+1.0}_{- 0.9}$
σcyc	keV	3.5 ± 0.6	$2.5^{+ 0.7}_{- 0.5}$	$3^{+4}_{-2}$	$2.8^{+1.0}_{- 0.5}$	$2.5^{+ 0.7}_{- 0.6}$	>2.12	$3.2^{+ 0.8}_{- 0.6}$	$0.8^{+1.0}_{- 0.6}$	<3.14
τcyc		$1.1^{+ 0.1}_{- 0.3}$	$0.9^{+ 0.2}_{- 0.3}$	$0.3^{+ 0.2}_{- 0.1}$	$0.9^{+ 0.2}_{- 0.4}$	$0.8^{+ 0.2}_{- 0.4}$	0.3 ± 0.1	$0.6^{+ 0.4}_{- 0.3}$	>0.59	$0.5^{+90.0}_{- 0.2}$
$\chi ^{2}_{\rm red}$(dof)		1.02 (434)	1.04 (422)	1.09 (423)	1.10 (422)	1.07 (425)	1.08 (421)	1.05 (425)	1.19 (413)	1.06 (419)
fdcut
Flux	× 10−10 erg cm−2 s−1	4.28 ± 0.03	5.48 ± 0.05	$4.57^{+ 0.03}_{- 0.05}$	4.19 ± 0.05	4.3 ± 0.2	8.7 ± 0.1	4.62 ± 0.05	$2.76^{+ 0.10}_{- 0.09}$	3.6 ± 0.1
NH	× 1022 cm−2	2.21 ± 0.03	2.14 ± 0.04	$2.22^{+ 0.06}_{- 0.07}$	$1.95^{+ 0.05}_{- 0.06}$	3.2 ± 0.2	2.7 ± 0.1	$2.05^{+ 0.07}_{- 0.08}$	2.5 ± 0.2	3.1 ± 0.2
NH, pcf	× 1022 cm−2	...	...	...	...	16 ± 3	...	...	$7^{+4}_{-3}$	$8^{+4}_{-2}$
fpcf		...	...	...	...	0.59 ± 0.03	...	...	0.4 ± 0.1	0.40 ± 0.09
Γ		1.15 ± 0.02	1.17 ± 0.02	$1.15^{+ 0.05}_{- 0.06}$	1.17 ± 0.04	$1.07^{+ 0.07}_{- 0.09}$	$0.85^{+ 0.06}_{- 0.07}$	$0.99^{+ 0.06}_{- 0.07}$	$1.32^{+ 0.08}_{- 0.10}$	$1.34^{+ 0.07}_{- 0.08}$
A	× 10−2 photons cm−2 s−1	4.4 ± 0.1	5.7 ± 0.2	5.0 ± 0.2	4.5 ± 0.2	$3.9^{+ 0.6}_{- 0.5}$	4.4 ± 0.1	$4.4^{+ 0.4}_{- 0.2}$	$3.7^{+ 0.6}_{- 0.5}$	$5.0^{+ 0.7}_{- 0.6}$
Ecut	keV	27 ± 2	28 ± 1	$24^{+6}_{-5}$	27 ± 3	$26^{+2}_{-3}$	$22^{+5}_{-4}$	$19^{+5}_{-6}$	$26^{+2}_{-3}$	28 ± 3
Efold	keV	5 ± 1	$4.0^{+1.0}_{- 0.9}$	7 ± 2	6 ± 2	5 ± 1	$7^{+1}_{-2}$	9 ± 2	5 ± 2	5 ± 2
EKα	keV	6.43 ± 0.01	6.43 ± 0.02	6.44 ± 0.02	6.40 ± 0.02	6.41 ± 0.01	6.39 ± 0.02	6.40 ± 0.02	6.42 ± 0.04	6.42 ± 0.03
IKα	× 10−4 photons cm−2 s−1	3.5 ± 0.3	4.4 ± 0.7	3.7 ± 0.6	3.3 ± 0.6	4.9 ± 0.5	7 ± 1	3.0 ± 0.6	$1.2^{+ 0.5}_{- 0.4}$	1.9 ± 0.5
IKβ	× 10−5 photons cm−2 s−1	4 ± 3	<8.76	<9.03	<7.58	9 ± 5	<14.20	<11.33	<3.69	8 ± 5
Ecyc	keV	$23.2^{+ 0.7}_{- 0.5}$	$23.2^{+ 0.8}_{- 0.6}$	$23^{+2}_{-1}$	$22.6^{+ 0.8}_{- 0.6}$	$22.3^{+ 0.6}_{- 0.5}$	$24^{+2}_{-1}$	$22.1^{+ 0.8}_{- 0.7}$	$22.0^{+ 0.8}_{- 0.6}$	23 ± 1
σcyc	keV	3.4 ± 0.6	2.7 ± 0.7	$3^{+2}_{-1}$	$2.4^{+1.0}_{- 0.8}$	$2.4^{+ 0.8}_{- 0.7}$	$3^{+2}_{-1}$	2 ± 1	<2.41	3 ± 1
τcyc		$0.9^{+ 0.2}_{- 0.1}$	0.8 ± 0.1	$0.4^{+ 0.2}_{- 0.1}$	$0.6^{+ 0.2}_{- 0.1}$	0.7 ± 0.1	$0.5^{+ 0.2}_{- 0.1}$	$0.5^{+ 0.6}_{- 0.1}$	>0.74	0.7 ± 0.2
$\chi ^{2}_{\rm red}$(dof)		1.04 (434)	1.06 (422)	1.11 (423)	1.10 (422)	1.09 (425)	1.08 (420)	1.03 (425)	1.20 (413)	1.07 (419)
npex
Flux	× 10−10 erg cm−2 s−1	4.24 ± 0.03	5.51 ± 0.06	4.51 ± 0.05	4.15 ± 0.05	$4.7^{+ 0.4}_{- 0.3}$	8.6 ± 0.1	4.55 ± 0.05	2.8 ± 0.2	$4.0^{+ 0.5}_{- 0.3}$
NH	× 1022 cm−2	2.05 ± 0.05	2.14 ± 0.07	$2.03^{+ 0.08}_{- 0.07}$	1.76 ± 0.07	3.4 ± 0.3	2.5 ± 0.1	1.83 ± 0.07	2.5 ± 0.2	3.2 ± 0.3
NH, pcf	× 1022 cm−2	...	...	...	...	19 ± 3	...	...	10 ± 4	$12^{+4}_{-3}$
fpcf		...	...	...	...	$0.66^{+ 0.05}_{- 0.06}$	...	...	$0.4^{+ 0.1}_{- 0.2}$	0.5 ± 0.1
α		$0.63^{+ 0.04}_{- 0.03}$	$0.79^{+ 0.06}_{- 0.05}$	$0.63^{+ 0.06}_{- 0.05}$	0.63 ± 0.05	0.8 ± 0.2	$0.36^{+ 0.07}_{- 0.06}$	$0.48^{+ 0.05}_{- 0.04}$	0.8 ± 0.2	1.1 ± 0.3
A	× 10−2 photons cm−2 s−1	4.3 ± 0.2	6.7 ± 0.5	$4.5^{+ 0.4}_{- 0.3}$	4.1 ± 0.3	$6^{+4}_{-2}$	4.3 ± 0.2	$3.5^{+ 0.3}_{- 0.2}$	$4^{+2}_{-1}$	$8^{+7}_{-3}$
B	× 10−3	3.2 ± 0.3	3.4 ± 0.3	2.4 ± 0.4	2.5 ± 0.4	2.8 ± 0.7	$4^{+1}_{-2}$	2.4 ± 0.4	$2.0^{+ 0.5}_{- 0.4}$	1.5 ± 0.5
kBT	keV	4.8 ± 0.1	4.7 ± 0.1	$5.2^{+ 0.4}_{- 0.2}$	5.1 ± 0.2	4.8 ± 0.2	$5.4^{+1.0}_{- 0.4}$	$5.5^{+ 0.3}_{- 0.2}$	4.7 ± 0.3	4.7 ± 0.3
EKα	keV	6.43 ± 0.01	6.43 ± 0.02	6.44 ± 0.02	6.40 ± 0.02	6.41 ± 0.01	6.39 ± 0.02	6.40 ± 0.02	6.42 ± 0.04	6.42 ± 0.03
IKα	× 10−4 photons cm−2 s−1	$3.6^{+ 0.4}_{- 0.3}$	4.6 ± 0.7	3.8 ± 0.6	3.4 ± 0.6	5.1 ± 0.6	8 ± 1	3.1 ± 0.6	1.2 ± 0.5	2.1 ± 0.6
IKβ	× 10−5 photons cm−2 s−1	4 ± 3	<8.80	<10.05	<8.45	10 ± 5	<14.67	<12.49	<4.18	$10^{+6}_{-5}$
Ecyc	keV	22.6 ± 0.5	22.3 ± 0.6	$23^{+2}_{-1}$	$22.3^{+ 0.7}_{- 0.6}$	21.9 ± 0.5	$25^{+4}_{-2}$	$22.3^{+ 0.8}_{- 0.7}$	$21.8^{+ 0.7}_{- 0.6}$	$22.6^{+1.0}_{- 0.9}$
σcyc	keV	2.7 ± 0.5	2.4 ± 0.7	$3^{+2}_{-1}$	$2.3^{+ 0.8}_{- 0.7}$	2.3 ± 0.7	$3^{+3}_{-2}$	2.0 ± 0.9	$1.0^{+1.0}_{- 0.7}$	3 ± 1
τcyc		$0.66^{+ 0.09}_{- 0.08}$	0.6 ± 0.1	0.4 ± 0.1	$0.6^{+ 0.2}_{- 0.1}$	0.6 ± 0.1	$0.5^{+ 0.4}_{- 0.1}$	$0.5^{+ 0.3}_{- 0.1}$	>0.69	$0.6^{+ 0.2}_{- 0.1}$
$\chi ^{2}_{\rm red}$(dof)		1.08 (434)	1.00 (422)	1.12 (423)	1.11 (422)	1.08 (425)	1.07 (420)	1.01 (425)	1.20 (413)	1.06 (419)

Notes. ^a2–10 keV unabsorbed flux. ^bPower-law normalization at 1 keV. ^cIron Kα intensity. ^dIron Kβ intensity.

Download table as:  ASCIITypeset image

This approximates the effect of the source being obscured by a varying absorbing column density. The entire source is absorbed by some base column density N_H (this corresponds to the value of N_H measured by the single tbnew component used in the other spectra), while a region with a higher column density of N_{H, pcf} covers some fraction f_pcf of the light from the source. This modification of the model is quite successful, bringing the value of reduced χ² for the dip spectrum and spectrum 9 down to ∼1 for each continuum model. The lowest-luminosity spectrum, number 8, still has the highest overall χ² even with the partial-covering model applied, with $\chi ^{2}_{\rm red} = 1.20$ for 413 degrees of freedom. We attempted to apply the partial-covering model to the remaining spectra, but the covering fraction could not be constrained and the quality of the fit was not improved. In the final two spectra, the addition of the partial-covering model significantly increases the power-law index relative to its value when only a single absorber is used, and the high-energy cutoff parameters are found at values that are more in line with their pre-dip values. While this means comparing the parameters of these spectra to those without the partial-covering model is somewhat risky, it should be noted that our pulse-to-pulse analysis (Section 4.3) finds broadly similar results in terms of N_H and the power-law index.

The fitted parameters for the three continuum models are listed in Table 2. We plot selected parameters against time and unabsorbed flux in Figure 4; as the different continuum models show similar behavior, and in the interest of clarity, we only plot the results from fdcut. For each plotted parameter versus flux, we perform a linear fit and calculate 90% error bars for the slope, as well as computing Pearson's product-moment correlation coefficient for the data; these results are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Linear Fits (with 90% Error Bars) and Correlation Coefficients for Time-resolved Spectroscopy

Parameter	Units	highecut			fdcut			npex
		Slopea	rb	p-valuec	Slope	r	p-value	Slope	r	p-value
Γ		$-0.054^{+ 0.006}_{- 0.006}$	−0.868	0.002	$-0.05^{+ 0.02}_{- 0.02}$	−0.846	0.004	...	...	...
Ecut	keV	$-1.3^{+ 0.4}_{- 0.4}$	−0.369	>0.1	$-0.5^{+ 0.9}_{- 0.9}$	−0.465	>0.1	...	...	...
Efold	keV	$0.6^{+ 0.3}_{- 0.3}$	0.078	>0.1	$0.2^{+ 0.3}_{- 0.3}$	0.327	>0.1	...	...	...
α		...	...	...	...	...	...	$-0.05^{+ 0.02}_{- 0.02}$	−0.689	0.040
kBT	keV	...	...	...	...	...	...	$0.06^{+ 0.06}_{- 0.06}$	0.599	0.088
Ecyc	keV	$0.3^{+ 0.2}_{- 0.2}$	0.586	0.097	$0.3^{+ 0.2}_{- 0.2}$	0.690	0.040	$0.4^{+ 0.2}_{- 0.2}$	0.896	0.001
σcyc	keV	$0.1^{+ 0.5}_{- 0.6}$	0.094	>0.1	$0.2^{+ 0.3}_{- 0.3}$	0.407	>0.1	$0.3^{+ 0.3}_{- 0.3}$	0.624	0.072
τcyc		$-0.12^{+ 0.04}_{- 0.04}$	−0.340	>0.1	$-0.03^{+ 0.04}_{- 0.04}$	−0.422	>0.1	$-0.02^{+ 0.03}_{- 0.03}$	−0.534	>0.1

Notes. ^aFrom linear fit against unabsorbed 2–10 keV flux in units of 10⁻¹⁰ erg cm⁻² s⁻¹. ^bPearson product-moment correlation coefficient. ^cTwo-tailed p-value for N = 9.

Download table as:  ASCIITypeset image

To obtain a more complete picture of the variability of the absorbing column and the iron line, we further divided the observation and extracted spectra for individual pulses. We only examine spectra for full pulses—spectra from fractional pulses are excluded to avoid contamination by the significant spectral variability within each pulse (see Section 4.4). After excluding fractional pulses, we have a total of 76 spectra. Only XIS spectra are analyzed, as the low exposure for each spectrum precludes the use of PIN data. Each spectrum is fit with a simple model consisting of a single absorbed power law with an additive Gaussian modeling the iron Kα line. This model generally fits the spectra quite well, with an average reduced χ² of 1.17. No evidence for partial covering was visible in the residuals, and the partial-covering model used in the previous section for the dip could not be constrained, so the absorption was modeled by a single tbnew component. The cflux model component in XSPEC was used to obtain unabsorbed fluxes for the power law and the iron line. We plot each parameter with respect to time and luminosity in Figure 5. We also calculated Pearson's correlation coefficients and performed linear fits on selected parameter versus the 2–10 keV power-law flux; these results are summarized in Table 4. The measured power-law index behaves roughly similar to what was seen in the time-resolved spectroscopy in Section 4.2, despite the use of a partial-covering model in that analysis. However, while the measured absorbing column density in the dip roughly lines up with the partial-covering N_H found in the time-resolved spectroscopy, the slightly elevated N_H in the final ∼10 pulse-to-pulse spectra do not reach the same levels as the partial-covering model found for the final two time-resolved spectra.

Zoom In Zoom Out Reset image size

The pulse profile of 4U 1538−522 is double-peaked, with a large primary pulse and a significantly smaller secondary pulse (see Figure 5), similar to that found by Robba et al. (2001), Rodes-Roca et al. (2010), and Hemphill et al. (2013). Our phase-resolved analysis was carried out by defining good time intervals for six phase bins, with phase zero defined as the peak of the main pulse. The remaining phase bins thus cover the rising and falling edges of the main pulse, the peak of the secondary pulse, and the low-flux dips on either side of the secondary pulse. This choice of phase binning covers all the major features of the pulse profile without sacrificing signal-to-noise. XIS and PIN spectra were analyzed simultaneously, using the three models described in Section 4.2. The parameters of these three models with 90% error bars are listed in Table 5, and the parameters for the fdcut continuum are plotted in Figure 6. The shape of the continuum, including the cyclotron line, changes considerably between the rising and falling edges of the main pulse. Interestingly, the absorbing column density appears to be significantly dependent on pulse phase (although the magnitude of the effect is relatively small, on the order of 10%). A physical mechanism for this effect is difficult to imagine, and on the whole unlikely—the absorber does not readily show signs of significant ionization, which one would expect if it was close enough to the neutron star to be affected by the star's rotation. It is likely that the dependence of N_H on pulse phase is in part due to our unphysical modeling of the continuum. This explanation is favored by the fact that the dependence of N_H on pulse phase varies with the choice of continuum model—under highecut, N_H is roughly constant in all bins except for the peak of the main pulse, where it drops suddenly, while in fdcut and npex the N_H follows a more "sawtooth" pattern, dropping sharply at the main pulse but rising slowly afterward. Meanwhile, the CRSF parameters and iron line intensity all show strong dependence on phase, without any strong dependence on continuum model.

Zoom In Zoom Out Reset image size

Table 5. Spectral Fits for Phase-resolved Spectroscopy

Parameter	Units	Phase Bin
		1	2	3	4	5	6
highecut
Fluxa	× 10−10 erg cm−2 s−1	6.32 ± 0.04	4.59 ± 0.04	2.98 ± 0.03	3.67 ± 0.03	2.60 ± 0.03	5.47 ± 0.04
NH	× 1022 cm−2	2.00 ± 0.03	2.34 ± 0.04	2.30 ± 0.05	2.35 ± 0.05	2.53 ± 0.05	2.35 ± 0.04
Γ		0.91 ± 0.01	1.13 ± 0.02	1.02 ± 0.02	0.98 ± 0.02	1.28 ± 0.02	1.18 ± 0.02
Ab	× 10−2 photons cm−2 s−1	4.2 ± 0.1	4.4 ± 0.1	2.37 ± 0.09	2.72 ± 0.09	3.2 ± 0.1	$5.7^{+ 0.2}_{- 0.1}$
Ecut	keV	23 ± 3	$12.5^{+ 0.6}_{- 0.7}$	$14^{+4}_{-1}$	$23^{+4}_{-8}$	$14^{+9}_{-2}$	$18^{+8}_{-4}$
Efold	keV	$6^{+1}_{-2}$	10.1 ± 0.7	$9.1^{+ 0.9}_{-2.0}$	$4^{+4}_{-1}$	$21^{+5}_{-10}$	$17^{+5}_{-7}$
EKα	keV	6.41 ± 0.02	6.41 ± 0.02	6.42 ± 0.02	6.43 ± 0.01	$6.415^{+ 0.009}_{- 0.008}$	6.42 ± 0.01
IKαc	× 10−4 photons cm−2 s−1	3.2 ± 0.5	3.3 ± 0.4	2.6 ± 0.4	4.0 ± 0.4	4.6 ± 0.4	3.9 ± 0.5
IKβd	× 10−5 photons cm−2 s−1	<8.00	<7.00	5 ± 4	5 ± 4	6 ± 3	<7.00
Ecyc	keV	$23.4^{+1.0}_{- 0.8}$	21.9 ± 0.6	21.9 ± 0.5	22 ± 2	21 ± 1	$21.4^{+ 0.7}_{- 0.6}$
σcyc	keV	$3.0^{+ 0.6}_{- 0.3}$	0.33 (fixed)	2 ± 1	$3.5^{+1.0}_{- 0.4}$	2.82 (fixed)	$3.0^{+ 0.6}_{-2.0}$
τcyc		0.9 ± 0.2	$2^{+2}_{-1}$	$1.0^{+2.0}_{- 0.3}$	$1.7^{+ 0.2}_{- 0.9}$	$0.3^{+ 0.4}_{- 0.1}$	$0.4^{+ 0.3}_{- 0.2}$
$\chi ^{2}_{\rm red}$(dof)		1.14 (432)	1.36 (423)	1.28 (417)	1.08 (423)	1.14 (427)	1.17 (438)
fdcut
Flux	× 10−10 erg cm−2 s−1	6.32 ± 0.04	$4.49^{+ 0.04}_{- 0.03}$	2.94 ± 0.03	3.63 ± 0.04	2.58 ± 0.03	5.46 ± 0.04
NH	× 1022 cm−2	1.99 ± 0.03	$2.00^{+ 0.06}_{- 0.04}$	2.17 ± 0.08	$2.24^{+ 0.07}_{- 0.08}$	$2.35^{+ 0.08}_{- 0.06}$	$2.31^{+ 0.04}_{- 0.05}$
Γ		$0.91^{+ 0.01}_{- 0.02}$	$0.72^{+ 0.08}_{- 0.03}$	$0.88^{+ 0.06}_{- 0.08}$	$0.86^{+ 0.06}_{- 0.07}$	$1.07^{+ 0.08}_{- 0.04}$	$1.14^{+ 0.03}_{- 0.04}$
A	× 10−2 photons cm−2 s−1	4.2 ± 0.1	$6.4^{+ 0.2}_{-1.0}$	$2.26^{+ 0.10}_{- 0.09}$	$2.56^{+ 0.10}_{- 0.09}$	$5.4^{+ 0.3}_{-2.0}$	5.7 ± 0.2
Ecut	keV	28 ± 2	<6.01	$16^{+3}_{-4}$	17 ± 3	<20.01	$29^{+3}_{-4}$
Efold	keV	3 ± 1	$8.8^{+ 0.2}_{- 0.5}$	$6.6^{+1.0}_{- 0.9}$	$6.0^{+ 0.7}_{- 0.9}$	$16.3^{+ 0.9}_{-2.0}$	8 ± 2
EKα	keV	6.41 ± 0.02	6.42 ± 0.02	6.42 ± 0.02	6.43 ± 0.01	6.415 ± 0.009	$6.42^{+ 0.02}_{- 0.01}$
IKα	× 10−4 photons cm−2 s−1	3.1 ± 0.5	$3.0^{+1.0}_{- 0.4}$	2.5 ± 0.4	3.9 ± 0.4	$4.5^{+3.0}_{- 0.4}$	3.8 ± 0.5
IKβ	× 10−5 photons cm−2 s−1	<8.00	<6.00	5 ± 4	<8.00	$6^{+5}_{-3}$	<7.00
Ecyc	keV	$24.1^{+ 0.8}_{- 0.5}$	$22.2^{+ 0.8}_{- 0.5}$	$22.2^{+ 0.5}_{- 0.4}$	$21.6^{+ 0.5}_{- 0.4}$	$21.9^{+1.0}_{- 0.9}$	$22.7^{+ 0.9}_{- 0.7}$
σcyc	keV	3.3 ± 0.6	0.16 (fixed)	1.5 ± 0.7	2.5 ± 0.6	2 ± 1	3 ± 1
τcyc		0.8 ± 0.2	>7.69	$1.2^{+2.0}_{- 0.3}$	$0.9^{+ 0.2}_{- 0.1}$	$0.4^{+ 0.8}_{- 0.1}$	$0.42^{+ 0.10}_{- 0.08}$
$\chi ^{2}_{\rm red}$(dof)		1.15 (432)	1.05 (423)	1.26 (417)	1.05 (423)	1.10 (426)	1.18 (438)
npex
Flux	× 10−10 erg cm−2 s−1	$6.31^{+ 0.05}_{- 0.04}$	4.42 ± 0.04	2.90 ± 0.03	3.60 ± 0.04	2.54 ± 0.03	5.40 ± 0.04
NH	× 1022 cm−2	1.91 ± 0.05	1.83 ± 0.05	1.95 ± 0.08	2.05 ± 0.08	$2.15^{+ 0.07}_{- 0.06}$	2.13 ± 0.05
α		0.50 ± 0.04	0.32 ± 0.03	0.33 ± 0.06	0.33 ± 0.06	0.64 ± 0.04	0.66 ± 0.03
A	× 10−2 photons cm−2 s−1	$4.5^{+ 0.3}_{- 0.2}$	$3.1^{+ 0.2}_{- 0.1}$	1.9 ± 0.2	2.4 ± 0.2	$2.4^{+ 0.2}_{- 0.1}$	5.1 ± 0.3
B	× 10−3	5.2 ± 0.3	2.8 ± 0.6	4.9 ± 0.8	6.7 ± 0.9	$1.1^{+ 0.4}_{- 0.3}$	1.8 ± 0.2
kBT	keV	5.0 ± 0.1	4.8 ± 0.2	4.6 ± 0.2	4.3 ± 0.2	$6.1^{+ 0.6}_{- 0.4}$	5.9 ± 0.2
EKα	keV	6.41 ± 0.02	6.42 ± 0.02	6.42 ± 0.02	6.43 ± 0.01	$6.415^{+ 0.009}_{- 0.008}$	6.42 ± 0.01
IKα	× 10−4 photons cm−2 s−1	3.2 ± 0.5	3.2 ± 0.5	2.6 ± 0.4	3.9 ± 0.5	4.6 ± 0.4	4.0 ± 0.5
IKβ	× 10−5 photons cm−2 s−1	<7.00	<8.00	5 ± 4	5 ± 4	7 ± 3	<8.00
Ecyc	keV	$23.9^{+ 0.5}_{- 0.4}$	$22.1^{+ 0.5}_{- 0.4}$	$22.3^{+ 0.5}_{- 0.4}$	$21.7^{+ 0.5}_{- 0.4}$	$21.9^{+1.0}_{- 0.8}$	22.4 ± 0.6
σcyc	keV	2.1 ± 0.5	0.33 (fixed)	$1.5^{+ 0.6}_{- 0.7}$	2.4 ± 0.5	2 ± 1	3.4 ± 0.7
τcyc		$0.48^{+ 0.09}_{- 0.08}$	$4^{+4}_{-2}$	$1.2^{+2.0}_{- 0.3}$	$0.9^{+ 0.2}_{- 0.1}$	$0.4^{+ 0.2}_{- 0.1}$	0.46 ± 0.06
$\chi ^{2}_{\rm red}$(dof)		1.24 (432)	1.05 (424)	1.27 (417)	1.07 (423)	1.07 (426)	1.15 (438)

Notes. ^a2–10 keV unabsorbed flux. ^bPower-law normalization at 1 keV. ^cIron Kα intensity. ^dIron Kβ intensity.

Download table as:  ASCIITypeset image

As can be seen in Figures 4 and 5, the column density and continuum parameters during the pre-dip portion of the observation stay relatively constant, with N_H at an average value of 2.16 × 10²² cm⁻² in the pulse-to-pulse data set. However, with the onset of the dip ∼32 ks into the observation, the measured column density increases drastically, first spiking to (24 ± 2) × 10²² cm⁻², quickly falling back to pre-dip levels, and then, immediately before the peak of the flare, rising again to (16 ± 2) × 10²² cm⁻². We note that the partial-covering model in the time-resolved spectrum of the dip reaches similar column densities ((19 ± 3) × 10²² cm⁻²) to the pulse-to-pulse results, suggesting that the necessity of the partial coverer may be due to the averaging of spectra with varying N_H, rather than a spatial partial coverer. The spikes in N_H suggest that the dip is caused by some overdense region of the accreted stellar wind passing through the line of sight. The duration of the increases in N_H along with the fact that the other continuum parameters do not change significantly during the dip indicates that the material producing the additional absorption must be relatively far from the star, as otherwise the material would be accreted onto the star in a relatively short time, changing the continuum emission and increasing the flux.

The continuum parameters (power-law index and cutoff and folding energies) are generally flat for the first half of the observation, remaining this way until a few pulsations before the peak of the flare, when the power-law index drops (in both the time-resolved spectra and the pulse-to-pulse spectra) by a factor of approximately 1/3. While the peak of the flare is clearly visible in the light curve, the point where the power-law index drops is a better approximation of the "true" starting point of the flare, with the beginning stages of the flare being obscured by the second spike in N_H. The unabsorbed flux is indeed climbing during the three pulses before the peak of the flare, as can be seen in the overplotted light curve in Figure 5. While the cutoff and folding energies do change slightly with the onset of the flare, they are poorly constrained compared to the power-law index, so it is difficult to see the effect of the flare in these parameters. The cutoff energy does appear to drop during the flare and in the first post-flare spectrum, before rising back to its pre-flare levels, while the even more poorly constrained folding energy (and kT in the npex model) rises during the flare. However, these effects are statistically very weak.

The source after the flare is less luminous, in terms of unabsorbed flux, than it was before the flare. In the pulse-to-pulse spectra, the power-law index additionally shows an overall negative correlation with luminosity, with a slope of −0.067 ± 0.008 per 10⁻¹⁰ erg cm⁻² s⁻¹. Pearson's coefficient for this is −0.58, which, for 76 points, corresponds to a two-tailed p-value of less than 10⁻⁵. This trend is preserved even if we exclude the dip and the spike from the calculation, although the slope drops slightly to −0.058 ± 0.011 per 10⁻¹⁰ erg cm⁻² s⁻¹. The time-resolved spectra show a similar correlation, with a slope of −0.05 ± 0.02 per 10⁻¹⁰ erg cm⁻² s⁻¹. This lends some credence to the partial-covering model used in the final two spectra of the time-resolved analysis—without the partial-covering model, the power-law index in the time-resolved spectra does not show a significant correlation with flux.

The continuum parameters generally show variability with phase, with the different continuum models displaying similar behavior. The power-law index, cutoff energy, folding energy, and npex's kT all show a peak somewhere around the rising edge of the main pulse, and roughly trend downward until the secondary pulse. The spectral shape of the rising and falling edges of the main pulse are distinctly different. Our continuum results are similar to those found using other satellites by, e.g., Clark et al. (1990), Coburn (2001), and Robba et al. (2001).

One potential area of concern is that the dip or flare represent some anomalous behavior of the source, and including those data with the non-flaring, less-obscured data may skew our results with regard to any observed correlations. However, the observed trends all persist even if we remove the dip and the flare from the sample, indicating that the trends are somewhat fundamental to the longer-term behavior of the source. The plot of iron line flux versus power-law flux with the dip and flare excluded has the same slope of (1.1 ± 0.2) × 10⁻³ as found before, while the power-law index has a slightly smaller slope, −0.058 ± 0.011, but still consistent with the earlier result. Removing the dip and flare obviously produces larger changes in the behavior of N_H, but the overall appearance of the parameter when plotted against the power-law flux is somewhat unchanged—there is an overall negative trend, as one would expect, with increased scatter in the measured N_H as the luminosity decreases. This increased scatter in N_H can be investigated by looking at the (lack of a) correlation between N_H and the power-law index: periods of increased N_H do not tend to be immediately related to changes in spectral hardness.

5.1.1. Connection between Absorption Event and Flare

The presence of the flare in the light curve is evidence of some degree of clumpiness or structure in the stellar wind of QV Nor (although simulations have shown that highly variable accretion can occur even in the presence of a relatively smooth stellar wind; see, e.g., Manousakis et al. 2014). The increased N_H seen before the flare is a possible source of the material accreted during the flare. However, if this is the case, the total amount of matter in the clump should reflect the observed luminosity during the flare. The N_H variability has a two-peak structure, with one large peak reaching ∼25 × 10²² cm⁻² several kiloseconds before the flare and a smaller peak immediately before the flare that reaches ∼15 × 10²² cm⁻². By following a method similar to that applied to Vela X–1 by Martínez-Núñez et al. (2014), we can estimate the size and, more importantly, the mass of the overdensity in the stellar wind. We will first assume that the spikes in N_H are caused by overdensities that move with the stellar wind. Using the orbital solution from Mukherjee et al. (2006), the Suzaku observation was carried out at an orbital phase of ∼0.3, and so we will additionally assume that the wind is moving roughly perpendicular to the line of sight.

In general, the wind velocity at a given distance will depend on the terminal wind velocity v_∞ via

$$\begin{equation} v_{\rm wind}(r) = v_{\infty } \left(1 - \frac{R_{\rm opt}}{r} \right)^{\beta }. \end{equation} \tag{ 7 }$$

The terminal wind velocity can in turn be estimated from the escape velocity for QV Nor via the correlation plot in Abbott (1982). From Clark et al. (1994), we have

$$\begin{eqnarray} v_{\rm esc} &=& [2GM_{\rm opt}(1-L_{\rm opt}/L_{\rm Edd})/R_{\rm opt}]^{1/2} \\ &\approx & 600\:\mathrm{km}\,\mathrm{s}^{-1}. \end{eqnarray} \tag{ 8 }$$

In these equations, M_opt ∼ 15 M_☉, R_opt ∼ 16 R_☉, β ∼ 0.9, L_opt = 1.6 × 10⁵ L_☉, and L_Edd ∼ 5 × 10⁵ L_☉ are the mass, radius, Roche lobe filling factor, luminosity, and Eddington luminosity of QV Nor (Reynolds et al. 1992; Rawls et al. 2011). From Abbott (1982), v_∞ ∼ 1300–1800 km s⁻¹. A range of values have been reported for the mass and radius of QV Nor, with Reynolds et al. (1992) finding M_opt = 19.9 ± 3.4 M_☉ and Rawls et al. (2011) finding 20.72 ± 2.27 M_☉ and 14.13 ± 2.78 M_☉ for the elliptical and circular orbital solutions, respectively, but these differences do not have large impacts on our results here due to the imprecision in estimating the terminal wind velocity.

At an inclination of ∼70 degrees (Rawls et al. 2011), the semimajor axis of the binary system is ∼1.7 × 10¹² cm, and at this distance the wind velocity is ∼7 × 10⁷ cm s⁻¹. The length scales of the overdensities can thus be estimated from their time of passage through the line of sight from v_windt_pass ∼ 2 × 10¹¹ cm for the large spike and ∼1 × 10¹¹ cm for the smaller spike. The average excess column densities over the extent of each spike in N_H are 9.1 × 10²² cm⁻² for the earlier, larger spike and 9.2 × 10²² cm⁻² for the second, smaller spike. The similarity in average column density is due to the shorter duration of the second spike (∼3 pulse periods) compared to the first spike (∼7 pulse periods). Given the estimated length scales, these imply number densities of ∼4 × 10¹¹ H atoms cm⁻³ and ∼8 × 10¹¹ H atoms cm⁻³, respectively, and thus masses of ∼10²² g and 2 × 10²¹ g. These represent a factor of ∼50–100 overdensity compared to the ∼5 × 10⁹ cm⁻³ in the wind at the neutron star's distance from QV Nor (Clark et al. 1994).

The X-ray flare persists for ∼7 pulse periods, if we define the "flare" state as when the source is spectrally harder (the green stars in Figure 5). If the material from the larger spike in N_H was fully accreted by the neutron star, the mass accretion rate would be ∼10¹⁸ g s⁻¹. If η, the efficiency of converting gravitational energy to X-rays, is ∼10%, the X-ray luminosity would be

$$\begin{eqnarray} L_{\rm flare} &=& \eta \frac{{ GM}_{\rm X}\dot{M}}{R_{\rm X}} \\ &\approx & 3 \times 10^{37}\:\mathrm{erg}\,\mathrm{s}^{-1}. \end{eqnarray} \tag{ 10 }$$

The 5–100 keV X-ray luminosity measured during the flare is somewhat lower than this, with a peak luminosity of 1 × 10³⁷ erg s⁻¹. However, the structure of the N_H variability, with the second peak during the flare and increased N_H post-flare, indicates that it is unlikely that the entirety of the overdense region was accreted onto the neutron star during the flare. Nonetheless, it is clear that the increased N_H seen during the dip represents enough material to support the observed flaring.

In both the time-resolved and the pulse-to-pulse spectra, the Fe Kα line has a significant positive correlation with flux. The points from the dip are slight outliers, remaining at roughly the same level as before the dip during the first spike in N_H, and rising with the flare during the second, smaller peak in N_H. The iron line's insensitivity to the dip is not surprising, considering that the flux above ∼5 keV is largely unaffected by the large increase in N_H that produces the dip. In the pulse-to-pulse spectra, Pearson's correlation coefficient is 0.701, which indicates a very strong correlation. A linear fit to the pulse-to-pulse data returns a slope of (1.1 ± 0.2) × 10⁻². The strong correlation between power-law flux and iron line flux implies that the iron must be fairly close to the neutron star (i.e., with a relatively small light travel time between the neutron star's continuum emission and the surrounding iron). The cross-correlation of the pulse-to-pulse iron line fluxes and power-law fluxes peaks at a lag of zero, which, at our time resolution, places an upper limit on this light travel time of 525.59 s, or a distance of ∼1.5 × 10¹³ cm. Meanwhile, phase-resolved spectroscopy shows the iron line varying with phase, following the general shape of the pulse profile, but with a significant phase lag of ∼0.67, or ∼350 s. If we assume that the iron emission peaks due to the light from the peak of the pulse, this would imply a distance of ∼10¹³ cm, or ∼1 AU. Work by Rodes-Roca et al. (2010) using XMM-Newton found variability of the iron line with phase, but they only report the iron line equivalent width, and with the complex continuum that they use, a direct comparison is difficult. Robba et al. (2001) divided the BeppoSAX spectrum of 4U 1538−522 into four phase bins, and while they could not constrain the intensity of the iron line especially well, their results are broadly in line with ours, with the brightest iron Kα offset from the main pulse.

Our constraints on the iron Kβ line provide some limited additional insight in this area, as the ratio of the Kβ intensity to the Kα intensity is a probe of the level of ionization of the iron. While we typically only find upper limits on the Kβ intensity, the values we do find are consistent with neutral iron, with the Kβ to Kα intensity ratio sitting around the expected 13.5% (Palmeri et al. 2005). This is again consistent with the distance found via the cross-correlation and the phase-resolved spectroscopy—if the iron is indeed at an average distance of ∼1 AU, we should not expect to see any strong signature of ionization.

The behavior of the cyclotron line over the course of the observation is our final feature of interest. This pseudo-absorption feature is produced when photons are scattered out of the line of sight by electrons in the accretion column above the neutron star's magnetic poles. In the ∼10¹² G magnetic field of the neutron star, the cyclotron motion of the electrons is quantized into Landau levels, and thus photons see an elevated scattering cross-section when they have sufficient energy to send an electron into the next level. This scattering results in a net reduction in line-of-sight photons at around the CRSF energy. The "12-B-12" rule gives an approximate relation between the observed CRSF energy and the magnetic field strength in the scattering region:

$$\begin{equation} E_{\mathrm{cyc}} = \frac{11.57\,\mathrm{keV}}{1+z} \times B_{12}, \end{equation} \tag{ 12 }$$

where B₁₂ is the magnetic field in units of 10¹² G and z ∼ 0.15 is the gravitational redshift in the scattering region. 4U 1538−522's CRSF energy of ∼23 keV implies a magnetic field strength of ∼2.2–2.3 × 10¹² G, assuming a radius of 10 km and a mass of either 0.87 ± 0.07 M_☉ or 1.104 ± 0.177 M_☉, per Rawls et al. (2011). Our measured cyclotron line energy is moderately higher than many previously measured values, with Clark et al. (1990), Robba et al. (2001), and Coburn (2001) all finding CRSF energies around ∼20–21 keV, while our CRSF energy is closer to ∼23 keV. The results of Rodes-Roca et al. (2009) and Hemphill et al. (2013) are closer in line with ours. While the possibility of a long-term change over time in the cyclotron line energy is intriguing, a proper comparison with past work would need to take much more careful account of the differing instruments and models (for both the continuum and the CRSF) used by different authors, as these may have not insignificant effects on the measured line energy (see, e.g., Müller et al. 2013).

The cyclotron line energy reaches a peak value of ∼24 keV during the flare, and there is a positive correlation between the CRSF energy and luminosity. Our linear fits to the data find a slope of 0.3 ± 0.2 keV per 10⁻¹⁰ erg cm⁻² s⁻¹ in the fdcut and highecut models and 0.4 ± 0.2 keV per 10⁻¹⁰ erg cm⁻² s⁻¹ when using npex. This is somewhat dependent on the inclusion of the point from the flare: if the flare is excluded, the slope for the fdcut continuum changes to 0.3 ± 0.3, which, while just barely consistent with zero, is still notably inconsistent with a negative slope.

Current theory regarding the cyclotron line, according to Becker et al. (2012), finds that the CRSF energy should exhibit different behavior with respect to source luminosity depending on the luminosity regime of the source. When the accretion column is locally super-Eddington (L > L_crit in Becker et al.), increases in luminosity should push the line-emitting region upward, toward lower magnetic fields, while below this critical luminosity, the opposite should occur. At still lower luminosities, the line-producing region should be relatively close to the neutron star surface, and is not expected to change drastically with changes in luminosity.

After computing the luminosity of 4U 1538−522 for each time-resolved spectrum from the 5–100 keV flux assuming a distance of 6.4 kpc (Reynolds et al. 1992), we can compare our results with the predictions of Becker et al. (2012). L_crit is given by Equation (32) in Becker et al. (2012):

$$\begin{eqnarray} L_{\rm crit} &=& 1.49 \times 10^{37}{\rm erg\,s}^{-1} \left(\frac{\Lambda }{0.1} \right)^{-7/5} w^{-28/15} \nonumber \\ &&\times \left(\frac{M}{1.4{ \,M}_{\odot }} \right)^{29/30} \left(\frac{R}{10{\rm \,km}} \right)^{1/10} \left(\frac{B_{\rm surf}}{10^{12}{\rm \,G}} \right)^{16/15}. \nonumber \\ \end{eqnarray} \tag{ 13 }$$

Here, M, R, and B are, respectively, the mass, radius, and surface magnetic field strength of the neutron star, w = 1 characterizes the shape of the photon spectrum inside the column (where there is assumed to be a mean photon energy of $\bar{E} = wkT_{\rm eff}$), and Λ characterizes the mode of accretion. Becker et al. only consider the case of Λ = 0.1, which approximates disk accretion. With this value for Λ, 4U 1538−522 is sub-critical, with the L_coul cut-off of Becker et al. (2012) sitting approximately between the highest non-flaring luminosity and the flare, and the Becker et al. (2012) predictions hold. However, 4U 1538−522 is believed to be a wind-accreting source, rather than a disk accretor, and thus we should choose Λ = 1 in our calculations. This presents us with an interesting result: the value for L_crit drops considerably, in fact, falling below the L_coul cut-off. It is uncertain as to what the theory predicts in this regime. This is the same as was found by Fürst et al. (2014) for Vela X−1, another wind-accreting HMXB of similar luminosity and with a CRSF of comparable energy (∼25 keV).

The problematic parameter here is Λ, which is difficult to determine precisely. The w parameter is unlikely to diverge much from w = 1, as the spectrum inside the accretion column is dominated by bremsstrahlung emission (Becker & Wolff 2007), and our results do not change significantly when w is varied around its assumed value of 1. The Λ parameter, on the other hand, is tied to the Alfvén radius R_A, as R_A in the disk of a disk-accreting source will be smaller than R_A for a wind-accretor experiencing a similar $\dot{M}$. Fürst et al. (2014) suggest that a very narrow accretion column could cause a breakdown of the theory in this case, and as the accretion column radius is tied to the Alfvén radius, differing values for Λ could have an effect here. However, it is also possible that the assumption of Λ = 1 is not precisely correct for 4U 1538−522; that is, the mode of accretion is not purely spherical. If the overdensity that produces the flare possesses significant angular momentum relative to the average angular momentum in the stellar wind, a transient disk may form during the flare and change the value of Λ from its value during the non-flaring state of the source.

To examine this, we fit the predicted E_cyc from Becker et al. (2012) to our cyclotron line energy measurements. The free parameters were Λ and the surface cyclotron line energy E_surf, although we limited E_surf to be greater than our largest measured E_cyc (i.e., we assumed that the magnetic field our observation samples is weaker than the magnetic field at the surface). The mass of the neutron star was set to 1.104 M_☉ (although using the lower-mass estimate from Rawls et al. 2011 did not significantly change our results), and the radius was assumed to be 10 km. As in Becker et al. (2012), we assume a dipole magnetic field. The theoretical E_cyc is thus dependent on the altitude h above the neutron star's surface:

$$\begin{equation} \frac{E_{\rm cyc}}{E_{\rm surf}} = \left(\frac{R + h}{R} \right)^{-3}, \end{equation} \tag{ 14 }$$

where h is given by

$$\begin{eqnarray} h_{\rm s} &=& 2.28 \times 10^{3}{\rm \,cm} \left(\frac{\xi }{0.01} \right)\left(\frac{M}{1.4{M}_{\odot }} \right)^{-1} \nonumber \\ && \times \left(\frac{R}{10{\rm \,km}} \right)\left(\frac{L_{\rm X}}{10^{37}{\rm \,erg\,s}^{-1}} \right) \end{eqnarray} \tag{ 15 }$$

for a supercritical source and by

$$\begin{eqnarray} h_{\rm c} &=& 1.48 \times 10^{5}{\rm \,cm} \left(\frac{\Lambda }{0.1} \right)^{-1} \left(\frac{\tau }{20} \right) \left(\frac{M}{1.4{M}_{\odot }} \right)^{19/14} \nonumber \\ && \times \left(\frac{R}{10{\rm \,km}}\right)^{1/14} \left(\frac{B_{\rm surf}}{10^{12}{\rm \,G}}\right)^{-4/7} \left(\frac{L_{\rm X}}{10^{37}{\rm \,erg\,s}^{-1}}\right)^{-5/7} \quad \end{eqnarray} \tag{ 16 }$$

for a subcritical source. This defines a piecewise function for E_cyc, which we fit to our results. For the fdcut continuum model, the fitted values were $\Lambda = 0.348^{+0.002}_{-0.058}$ and $E_{\rm surf} = 24.68^{+0.36}_{-0.03}$ keV. This best-fit result is plotted over the measured values for E_cyc in Figure 7. We additionally fit the predicted cyclotron line energy fixed Λ = 1.0 and Λ = 0.1. While the quality of the Λ = 0.348 fit is poor (χ² = 12.3 for 7 degrees of freedom), it is clearly superior to the pure-disk and pure-wind assumptions. The fitted function additionally implies a critical luminosity of ∼6 × 10³⁶ erg s⁻¹, in between the flaring and non-flaring states of 4U 1538−522. This would suggest that during the flare, the luminosity reached a sufficiently high level to create a stationary radiative shock in the accretion column, which pushes down the average height of the line-producing region. However, we lack data covering the transition from sub-critical to super-critical accretion, and the discontinuous change in this formulation is difficult to physically justify.

Zoom In Zoom Out Reset image size

Alternative treatments of the CRSF-luminosity correlation can be found in Poutanen et al. (2013) and Nishimura (2014). Poutanen et al. posit that the CRSF is formed out of light from the accretion column reflecting off of the neutron star surface, finding that, for a dipole field, a higher-altitude emission region in the accretion column will illuminate a larger fraction of the stellar surface and sample a weaker average magnetic field. This results in smaller predicted correlations compared to Becker et al. (2012), but the signs of those correlations are preserved, as the cyclotron line energy is still dependent on the altitude of some scattering region in the accretion column. Their work focuses on much higher luminosities than we see in 4U 1538−522, even during the flare, but their results should be somewhat compatible with ours, as their observed trends are still dependent on variations in the altitude of the emitting region in the accretion column. One area of concern is that the Poutanen et al. model suggests that the CRSF width should, to some extent, vary inversely with the CRSF energy, as lower-energy lines are produced by sampling a larger range of magnetic fields in the neutron star atmosphere, smearing out the line. However, this interpretation is simplistic, and additional study is needed to determine exactly what the reflection model predicts for the variation of line width with luminosity.

Meanwhile, Nishimura simulates a superposition of cyclotron lines produced in the accretion column within a range of altitudes and manages to reproduce the observed correlations in multiple sources. It is interesting to note here that, in the flare, the primary pulse is considerably brighter than the secondary, and in our phase-resolved analysis, the cyclotron line energy peaks with the primary pulse. This may be indicative of behavior similar to how Nishimura explains the observed positive correlation in Her X−1: the increased primary pulse relative to the rest of the pulse profile weighs that phase bin's CRSF energy above the others and produces the observed positive trend. However, this would imply that one should see a correlation between the ratio of the height of the primary pulse to the height of the secondary pulse. By modeling each double-peaked pulse with a sum of two Gaussians and computing the ratio of the heights, we can confirm that this is not the case—the average pulse height ratio for each time bin used in the time-resolved spectroscopy is not significantly correlated with E_cyc (see Figure 8).

Zoom In Zoom Out Reset image size

The CRSF width is also positively correlated with luminosity; this is in line with the observed correlation between CRSF energy and width seen in multiple sources and, indeed, across multiple sources (see, e.g., Coburn 2001)). The correlation between E_cyc and σ_cyc is linear, with a slope of 0.9 ± 0.6. The ratio of σ_cyc to E_cyc is constant at ∼0.11, with the exception of the lowest-luminosity measurement from time bin 8. The ratio of E_cyc to σ_cyc in a self-emitting atmosphere is given by Meszaros & Nagel (1985) as

$$\begin{equation} \frac{\sigma _{\rm cyc}}{E_{\rm cyc}} \approx \left(8 \ln 2 \frac{kT_{\rm e}}{m_{\rm e}c^{2}}\right)^{\frac{1}{2}}|\cos \theta |. \end{equation} \tag{ 17 }$$

As θ is the angle between the magnetic field and the line of sight, cos θ should not vary in phase-averaged spectra from a single observation. Thus, the interpretation is that a constant E_cyc/σ_cyc ratio implies a relatively constant electron temperature of ∼1.3 keV. The largest outlier is time bin 8, where the width of the CRSF drops to ∼1 keV. The line is easily detected in this spectrum, appearing clearly in the residuals, but its shape is not fit well by the gauabs model. The cyclabs model, which uses a Lorentzian line profile, likewise does not fit the line well. This stems from the fact that there is no physics-based model currently available for the CRSF, and so the fitted width in these bins may not be physically meaningful. The CRSF depth is relatively stable throughout the observation, within errors, although the low width in time bin 8 results in only a lower bound for the optical depth.

The cyclotron line varies significantly with phase (Figure 6), and remains in phase with the pulse profile in all continuum models. The behavior of the cyclotron line energy and width are comparable to the results of Clark et al. (1990) and Coburn (2001). The peak value of E_cyc is ∼24 keV, similar to the peak cyclotron line energy seen in the time-resolved spectra. The CRSF width is also correlated with the pulse profile, although this is a more conditional statement than it was for the CRSF energy. The falling edge of the primary pulse is particularly problematic, as the CRSF width in that phase bin is exceptionally low, at ∼0.1 − 0.3 keV. The depth in this bin is unbounded at the upper end due to the low width, and the line is overall much more poorly detected in comparison to the other phase bins, with its addition lowering the value of χ²/dof from 474.68/425 to 445.04/423 (compare to the peak of the main pulse, where the CRSF lowers χ²/dof from 689.74/435 to 495.89/432). However, beyond this one phase bin, the width shows a slight rise with phase, peaking with the primary pulse, and there is a hint of a correlation between the cyclotron line width and energy, similar to that seen in the phase-averaged spectra, but weaker. The depth follows the opposite pattern, with the highest-energy CRSF being one of the shallower features. The CRSF depth also displays the same asymmetry across the main pulse as seen in the other spectral parameters. Although Equation (17) predicts variation of σ_cyc/E_cyc with luminosity within the pulse, no significant correlation can be found due to the large error bars.