An XMM-Newton and NuSTAR Study of IGR J18214-1318: A Non-pulsating High-mass X-Ray Binary with a Neutron Star

NuSTAR is the first hard X-ray focusing telescope in space, providing 58'' half-power diameter angular resolution (Harrison et al. 2013). The NuSTAR FPMs cover the 3–79 keV band with moderate energy resolution (0.4 keV energy at 6 keV) and operate with 0.1 ms temporal resolution (Harrison et al. 2013). We processed the data from the two NuSTAR instruments, FPMA and FPMB, with the NuSTARDAS pipeline software v1.4.1, the 20150612 version of the NuSTAR Calibration Database (CALDB), and High Energy Astrophysics Software (HEASOFT) version 6.16.

Cleaned event lists were produced with the routine nupipeline. We used nuproducts to extract spectra, including response matrix files (RMFs) and ancillary response files (ARFs), and light curves, applying barycenter and dead-time corrections. In order to choose the aperture region sizes, we measured the surface brightness of the profile of the source by measuring the average count rate per pixel in concentric annuli centered on the source. We found that the profile substantially flattens at a distance of $150^{\prime\prime}$ from the source, so to guarantee a clean background measurement, we defined the background region as an annulus centered on the source with an inner radius of $200^{\prime\prime}$ and an outer radius of $250^{\prime\prime}$. Within $60^{\prime\prime}$ of the source, the count rate per pixel is at least 10 times higher than in the background region, so we defined the source aperture region as a circle with a $60^{\prime\prime}$ radius (corresponding to the PSF encircled energy fraction of 75%) to limit background contamination while still ensuring good photon statistics (∼30,000 counts in FPMA and FPMB combined). Figure 1 displays the NuSTAR FPMA observation with the source and background regions. We checked the light curves from the background regions for significant count rate variations, but did not find any significant background variability in the 3–79, 3–12, or 12–30 keV bands. The 3–30 keV dead-time-corrected source count rate is 1.2 counts s⁻¹ for FPMA and FPMB combined.

Zoom In Zoom Out Reset image size

We made use of data from the EPIC pn and MOS instruments on XMM-Newton, which provide $15^{\prime\prime}$ HPD angular resolution at soft X-ray energies. Observations were performed in small window mode with a medium filter. The EPIC pn instrument covers the 0.3–12 keV band with an energy resolution of 150 eV at 6 keV and, in the small window mode, it has 5.7 ms time resolution (Strüder et al. 2001). The MOS cameras provide a similar energy resolution in the 0.3–10 keV band but have a poorer timing resolution of 0.3 s in the small window mode (Turner et al. 2001).

We processed the XMM-Newton data with Science Analysis Software v13.5.0, making images, spectra, and light curves for EPIC pn, MOS1, and MOS2. To look for contamination from proton flares, we made EPIC pn and MOS light curves in the 10–12 keV bandpass, but we did not find any significant flares. For all three instruments, we made new event lists using the standard filtering criteria,¹² and we converted the photon arrival times to the solar system barycenter.

We used xmmselect to extract the data products, using the images primarily to create the source and background extraction regions. The source extraction region is a circle with a 40'' radius centered on the source position. For EPIC pn, the background region is a rectangle with an area of 1.4 square arcminutes that is located approximately 2' from the source. For the MOS detectors, rectangular background regions are also used, but they are farther away from the source because they had to be located on one of the outer MOS CCDs due to the CCD configuration of the small window mode. The 0.3–12 keV live-time-corrected count rate is 1.3 counts s⁻¹ with EPIC pn and 0.8 counts s⁻¹ for the two MOS cameras combined. We checked the observations for pile-up and found that it was not an issue.

First, we describe the spectral analysis of the flux-averaged source spectrum. We extracted spectra, ARFs, and RMFs from the XMM-Newton EPIC pn, MOS1/2, and NuSTAR FPMA/B instruments as described in Section 2. The spectra were rebinned with the requirement that the signal significance in each bin be $\geqslant 10$, except for the highest energy bin which was required to have a significance $\geqslant 3$. We used the XSPEC version 12.8.2 software to jointly fit the five spectra, allowing for different calibration constants for each instrument. The cross-calibration constants for the MOS1 and MOS2 instruments were consistent in all the fits, differing by less than 1σ from each other, so we linked the MOS1/2 constants together, removing one free parameter from the models.

We first fit the data using an absorbed power-law model (tbabs∗powerlaw), adopting the abundances from Wilms et al. (2000) and photoionization cross-sections from Verner et al. (1996). This simple model, which was sufficient for describing previously available soft X-ray data with lower photon statistics (Tomsick et al. 2008; Rodriguez et al. 2009), yields a poor fit (${\chi }_{\nu }^{2}=3.3$ for 585 dof). As can be seen in the residuals in Figure 8, this simple power-law fit underestimates the flux below 2 keV and overestimates it above 20 keV. Accounting for the flux above 20 keV requires introducing an exponential cutoff to the power-law spectrum, while the soft excess can be accounted for either by adding a blackbody component (Model 1) or a partial-covering absorber (Model 2), which provide equally good fits. Adding only one of these components (highecut, bbody, or pcfabs) to the absorbed power-law model is insufficient, leaving large residuals either below 2 keV or above 20 keV.

Zoom In Zoom Out Reset image size

The spectral fits and residuals resulting from our best-fitting models are shown in Figure 9. These models also include a Gaussian line to fit the Fe Kα line emission at 6.40 ± 0.02 keV, which is clearly visible in Figure 10. The energy of this line indicates it must originate in cool, low-ionization material located in the supergiant wind (Torrejón et al. 2010). The spectral parameters of the best-fitting models are listed in Table 2. As can be seen, the reduced ${\chi }^{2}$ values of the spectral fits are good enough that no additional components are required, and no prominent features remain in the residuals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Best-fit Spectral Parameters

	Model 1	Model 2
	tbabs∗(bbody+	tbabs∗pcfabs∗
	powerlaw∗highecut	(powerlaw∗highecut
	+Gaussian)	+Gaussian)
${N}_{{\rm{H}}}$ (1022 cm−2)	${4.2}_{-0.2}^{+0.3}$	4.3 ± 0.6
${{kT}}_{\mathrm{BB}}/{N}_{{\rm{H}},\mathrm{partial}}$	${1.74}_{-0.05}^{+0.04}$ keV	${9.8}_{-1.1}^{+1.5}\times {10}^{22}$ cm−2
BB norm./Cov. frac.	$1.3\pm 0.1\times {10}^{-11}$	${0.77}_{-0.06}^{+0.05}$
Γ	${0.4}_{-0.4}^{+0.3}$	${1.48}_{-0.07}^{+0.08}$
PL norm.	$4.9\pm 2.0\times {10}^{-12}$	$2.6\pm 0.2\times {10}^{-11}$
${E}_{\mathrm{cut}}$ (keV)	${12.0}_{-1.3}^{+1.0}$	${7.4}_{-0.5}^{+0.6}$
${E}_{\mathrm{fold}}$ (keV)	${14.0}_{-1.5}^{+3.2}$	${23.0}_{-2.4}^{+3.3}$
${E}_{\mathrm{line}}$ (keV)	${6.40}_{-0.02}^{+0.03}$	6.40 ± 0.02
${\sigma }_{\mathrm{line}}$ (eV)	<85	<102
EWline (eV)	${53}_{-21}^{+16}$	${57}_{-14}^{+16}$
${C}_{\mathrm{MOS}\mathrm{1,2}}$	1.04 ± 0.02	1.04 ± 0.02
${C}_{\mathrm{FPMA}}$	1.25 ± 0.02	1.24 ± 0.02
${C}_{\mathrm{FPMB}}$	1.32 ± 0.03	1.31 ± 0.03
${\chi }_{\nu }^{2}$/dof	1.10/578	1.10/578

Note. Errors provided are 90% confidence. Cross-normalizations between instruments are calculated relative to the XMM-Newton EPIC pn instrument. Abbreviations: BB—blackbody, PL—power law. The BB normalization is the unabsorbed 0.3–10 keV flux (in units of erg cm⁻² s⁻¹) of the BB component. The PL normalization is the 0.3–10 keV flux (in units of erg cm⁻² s⁻¹) of the PL component (in the case of the pcfabs model, it is the flux of the PL component that is subject to local absorption by the partial-covering absorber).

Download table as:  ASCIITypeset image

However, in order to statistically test for the presence of cyclotron lines, we added a cyclotron absorption component (cyclabs) to our models and performed new spectral fits. Since the cyclotron line width was very poorly constrained when left as a free parameter, we set its upper limit to 10 keV, since cyclotron line widths of accreting X-ray pulsars typically fall in the 1–10 keV range (Coburn et al. 2002). We also set the optical depth of the second harmonic to zero since it could not be constrained. The spectral parameters of all other Model 1 and 2 components were allowed to vary in order to find the best-fitting model that includes cyclabs as a multiplicative component. The best-fitting cyclotron line parameters derived by adding cyclabs to Model 1 are an optical depth ${\tau }_{\mathrm{cyc}}={0.25}_{-0.18}^{+0.21}$ and line energy ${E}_{\mathrm{cyc}}={27}_{-6}^{+4}$ keV; the inclusion of cyclabs only reduced the chi-squared value of the fit by 5 and left ${\chi }_{\nu }^{2}$ unchanged. The cyclotron parameters derived by adding cyclabs to Model 2 are ${\tau }_{\mathrm{cyc}}=0.16\pm 0.06$ and ${E}_{\mathrm{cyc}}={11}_{-7}^{+4}$ keV, which differs from the cyclotron energy found for Model 1; adding the cyclabs component to Model 2 improved the chi-squared value by 19.2 and reduced ${\chi }_{\nu }^{2}$ to 1.07 from 1.10, a marginal improvement on the quality of the fit.

In order to determine the significance of this improvement to the chi-squared value for Model 2, we generated 1000 simulated data sets, including both the NuSTAR and XMM data and followed the procedure applied in Bellm et al. (2014), Bhalerao et al. (2015), and Bodaghee et al. (2016). Each simulated data set was fit by the null model (Model 2 without cyclabs) and the test model with a cyclabs feature, and the difference in chi-squared values (${\rm{\Delta }}{\chi }^{2}$) between the two model fits was calculated. The maximum value of ${\rm{\Delta }}{\chi }^{2}$ from these simulations was 19.3, slightly higher than the observed value. Based on the distribution of ${\rm{\Delta }}{\chi }^{2}$ from our simulations, we estimate there is roughly a 0.001% chance of measuring the observed value of ${\rm{\Delta }}{\chi }^{2}=19.2$ by chance, and that therefore the significance of the cyclotron line in IGR J18214-1318 is about $3.3\sigma$. Given the fact that this detection is marginal and dependent on adopting Model 2 rather than Model 1 for the soft excess, it does not constitute substantive evidence for the presence of a cyclotron absorption feature.

Nonetheless, the absence of such features does not disprove the possibility that IGR J18214-1318 harbors an NS. The fact that the e-folding energy of the exponential cutoff is $\lt 25\,\mathrm{keV}$, regardless of which of the two best models is adopted, strongly suggests that IGR J18214-1318 is an NS HMXB. Furthermore, the photon index below the cutoff is harder when adopting the blackbody rather than the partial-covering model, but in both cases is within the range observed in NS HMXBs, which tend to exhibit harder photon indices than BH HMXBs (Coburn et al. 2002).

For Model 1, the blackbody component accounting for the soft excess has a temperature of 1.7 keV, which is higher than the ${kT}\approx 0.1\,\mathrm{keV}$ thermal component exhibited by BH HMXBs in the hard state (Di Salvo et al. 2001; McClintock & Remillard 2006; Makishima et al. 2008); BH HMXBs can exhibit blackbody temperatures as high as 2 keV in the soft state, but the power-law component of IGR J18214-1318 is much stronger than that of a BH in the soft state (McClintock & Remillard 2006). Assuming that IGR J18214-1318 lies at a distance of 9–10 kpc, as favored by the properties of its near-IR counterpart (Butler et al. 2009), the radius of the blackbody-emitting region is 0.3 km, which is consistent with the size of NS hot spots. However, while the blackbody interpretation thus provides some additional evidence in favor of the NS hypothesis, the soft excess seen in the spectrum can be equally well-fit by a partial-covering model. Using this model, we measure that the whole system lies behind a column density of 4 × 10²² cm⁻², which is just in excess of the Galactic interstellar column density integrated along the line of sight, and that about 77% of the X-ray emission is obscured by an additional column density of ∼10²³ cm⁻². This partial-covering absorber can be attributed to dense clumps in the supergiant wind, and thus does not provide any additional insight into the compact object in this HMXB.

The mean 3–12 keV flux is ${1.70}_{-0.05}^{+0.02}\times {10}^{-11}$ erg cm⁻² s⁻¹ in the XMM-Newton EPIC pn observations and ${2.18}_{-0.10}^{+0.03}\times {10}^{-11}$ erg cm⁻² s⁻¹ in the NuSTAR observations (FPMA and FPMB averaged). At a distance of 9–10 kpc, these fluxes correspond to unabsorbed luminosities of $\approx 1\mbox{--}2\times {10}^{35}$ erg s⁻¹. Through simultaneous NuSTAR and XMM-Newton observations of PKS 2155-304 and 3C 273, Madsen et al. (2015) found that the flux cross-calibrations of the XMM-Newton and NuSTAR instruments are accurate to better than 7%; in our spectral fits, the NuSTAR fluxes (see Table 2) are higher than the XMM-Newton fluxes by about 30%, primarily because the NuSTAR observations cover a longer duration of time and the source undergoes some large flares after the XMM-Newton observations end, as shown in Figure 6. The average X-ray flux during these observations is a factor of 3 lower than the average flux during the 2008 and 2009 observations (Tomsick et al. 2008; Rodriguez et al. 2009), and the light curves in Figure 6 show that on hour-long timescales, the flux can vary by more than a factor of 20.

In order to verify whether the flux-dependent spectral variations are small in IGR J18214-1318, as suggested by the stability of the hardness ratios with time, and to investigate the origin of the spectral variations if they exist, we made low flux state and high flux state spectra. The low (high) flux XMM-Newton EPIC pn and NuSTAR FPMA/B spectra consist of counts extracted from 100 s time intervals of simultaneous XMM-Newton and NuSTAR observations during which the NuSTAR, FPMA and FPMB combined, 3–12 keV count rate is less (greater) than 0.5 counts s⁻¹. This count rate threshold was chosen taking into account the point at which the hardness ratio changes most substantially in Figure 7 and trying to ensure that the low and high state spectra would have roughly the same number of total source counts. The spectra were rebinned so that the signal significance in each energy bin was $\geqslant 7$, except for the highest energy bin which was only required to have a significance $\geqslant 3$.

First, we fit the low and high state spectra independently with a simple absorbed power-law model. Both the low and high state spectra exhibit data residuals to the model fit that are very similar to those of the flux-averaged spectrum shown in Figure 8, with an excess below 2 keV and a deficit above 10 keV. In both cases, in order to remove the structure in these residuals and produce a good ${\chi }_{\nu }^{2}$ value, it is necessary to both add a high-energy cutoff to the power law and to introduce either a blackbody or a partial-covering absorber to model the soft excess. Thus, both the soft excess and high-energy cutoff appear to be persistent features of the spectrum of IGR J18214-1318. The addition of a Gaussian component does not significantly improve the ${\chi }_{\nu }^{2}$ value of the fits to either the low or high state spectra, as the sensitivity of the low and high state spectra is lower than that of the flux-averaged spectrum and insufficient to detect the relatively weak Fe line at 6.4 keV.

In order to establish whether any statistically significant spectral variations exist between the low and high state spectrum, we fit the low and high state spectra together, allowing only the cross-normalization for each spectrum to vary independently. Adopting Model 1 (constant∗tbabs∗(bbody+powerlaw∗highecut)) results in ${\chi }_{\nu }^{2}=1.05$ for 483 dof, while adopting Model 2 (constant∗tbabs∗pcfabs∗powerlaw∗highecut) results in ${\chi }_{\nu }^{2}=1.07$ for 483 dof. The probability that the data for both spectra are drawn from the same model is 22% for Model 1 and 14% for Model 2. Thus, the ${\chi }_{\nu }^{2}$ values are statistically consistent with a flux-independent spectral model.

However, looking in detail at the χ residuals to the flux-independent fits shown in Figure 11, it is striking that regardless of whether Model 1 or 2 is adopted, the lowest energy bins are offset by roughly 3σ from the model fit, with the models underestimating the high state spectrum and overestimating the low state spectrum below 2 keV. To assess whether this deviation at soft energies is significant, we rebinned the high state spectrum to match the energy binning of the low state spectrum so as to be able to compare the χ residuals of the two spectra on a bin-by-bin basis. We refit the high and low state spectra with Model 1 and Model 2, only allowing the cross-normalization for each spectrum to vary independently. The largest difference between the χ residuals of the high state and low state spectra occurs in the lowest energy bin. For Model 1(2), the lowest energy bin of the high state spectrum is 3.98σ (3.32σ) above the model prediction while the lowest energy bin of the low state spectrum is 2.72σ (3.31σ) below the model prediction, resulting in a net difference of 6.70σ (6.63σ). If a model is a good fit to the data, we expect the χ residuals of each spectrum to follow a Gaussian distribution with a standard deviation of 1.0; since the distribution of the difference of two normally distributed variables is also a Gaussian with variance equal to the sum of the variances of the two variables, for a good model fit we expect the difference of the χ residuals of the high state and low state spectra to follow a Gaussian distribution with a standard deviation of $\sqrt{2}$. Given that the spectra have 223 energy bins, the difference between the χ residuals of the two spectra should exceed 4.0σ in at most one bin. As shown by the pink histograms in Figure 12, regardless of whether Model 1 or 2 is used, when the spectral parameters are not allowed to vary between the high and low state, three energy bins have χ residual differences exceeding 4.0σ, and the difference of 6.6σ–6.7σ exhibited by the lowest energy bin is much higher than any other bin.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Therefore, we explore whether decoupling one or more spectral parameters between the high state and low state spectral models significantly improves the agreement between the observed and expected distributions of the χ residual differences. We find that for both Model 1 and Model 2 a minimum of two spectral parameters must be decoupled between the high state and low state spectral models in order for the χ residual difference to be $\lt 4\sigma$ in all but 1 of the 223 energy bins. For Model 1, the two parameters whose decoupling leads to the greatest improvement in ${\chi }_{\nu }^{2}$ are the hydrogen column density and the ratio of the blackbody flux to the power-law flux, while for Model 2 they are the covering fraction and the column density of the partial-covering absorber. The blue striped histograms in Figure 12 show the resulting χ residual difference distribution when these parameters are allowed to vary with source flux in the spectral models.

Table 3 presents the best-fit results for the flux-dependent models which allow the aforementioned parameters to vary between the high and low flux state. When adopting Model 1, the high state spectrum is better fit by a lower absorbing column density and a lower ratio of the blackbody flux to the power-law flux than the low state spectrum. The first trend may suggest that a higher fraction of the local absorbing medium is ionized as the X-ray flux increases, and the second may suggest that the size of the hot spot is larger than the accretion column, resulting in the blackbody emission not being as sensitive to fluctuations in the accretion rate as the power-law emission (Comptonized bremsstrahlung; Becker & Wolff 2007; Farinelli et al. 2016) from the accretion column. When adopting Model 2, the high state spectrum is better fit by a partial-covering absorber with a higher column density and a lower covering fraction than the low state spectrum. This trend is consistent with the accretion rate increasing as the compact object travels through regions where the stellar wind is clumpier.

Table 3.  Variations of Spectral Parameters with Source Brightness

	Model 1			Model 2
	tbabs∗(bbody+powerlaw∗highecut)			tbabs∗pcfabs∗powerlaw∗highecut
	Both	Low	High	Both	Low	High
${N}_{{\rm{H}}}$ (1022 cm−2)	⋯	${4.4}_{-0.4}^{+0.6}$	${3.6}_{-0.4}^{+0.7}$	4.0 ± 1.1	⋯	⋯
${{kT}}_{\mathrm{BB}}$ (keV)/${N}_{{\rm{H}},\mathrm{partial}}$ (1022 cm−2)	${1.74}_{-0.18}^{+0.06}$	⋯	⋯	⋯	7 ± 2	${11}_{-3}^{+4}$
BB norm. (10−11 erg cm−2 s−1)/Cov. frac.	⋯	1.07 ± 0.14	0.85 ± 0.15	⋯	0.9 ± 0.1	0.8 ± 0.1
Γ	${0.2}_{-0.6}^{+0.5}$	⋯	⋯	${1.51}_{-0.14}^{+0.10}$	⋯	⋯
PL norm. (10−11 erg cm−2 s−1)	${0.35}_{-0.13}^{+0.20}$	⋯	⋯	${2.3}_{-0.3}^{+0.2}$	⋯	⋯
${E}_{\mathrm{cut}}$ (keV)	${12}_{-4}^{+1}$	⋯	⋯	${7.3}_{-0.9}^{+0.7}$	⋯	⋯
${E}_{\mathrm{fold}}$ (keV)	${13}_{-3}^{+5}$	⋯	⋯	${23}_{-4}^{+5}$	⋯	⋯
${C}_{\mathrm{XMM},\mathrm{pn}}$	⋯	1.0 (fixed)	${2.68}_{-0.22}^{+0.24}$	⋯	1.0 (fixed)	${2.44}_{-0.11}^{+0.12}$
${C}_{\mathrm{FPMA}}$	⋯	1.07 ± 0.04	${2.80}_{-0.21}^{+0.24}$	⋯	1.10 ± 0.04	${2.63}_{-0.11}^{+0.12}$
${C}_{\mathrm{FPMB}}$	⋯	1.10 ± 0.05	${2.99}_{-0.23}^{+0.25}$	⋯	1.14 ± 0.05	${2.81}_{-0.12}^{+0.13}$
${\chi }_{\nu }^{2}$/dof	1.01/481	⋯	⋯	1.02/481	⋯

Note. The low (high) flux state includes all 100 s time intervals of simultaneous XMM and NuSTAR observations during which the 3–12 keV NuSTAR count rate is lower (higher) than 0.5 counts s⁻¹ (corresponding to a 3–12 keV XMM-Newton EPIC pn count rate of 1.3 counts s⁻¹). Errors provided are 90% confidence. Cross-normalizations between instruments are calculated relative to the XMM-Newton EPIC pn instrument low flux state spectrum. For each spectral model, the low and high flux state spectra are fit together, with the parameters in the "Both" column tied to one another. Abbreviations: BB—blackbody, PL—power law. The BB normalization is the unabsorbed 0.3–10 keV flux of the BB component. The PL normalization is the 0.3–10 keV flux of the PL component (in the case of the pcfabs model, it is the flux of the PL component that is subject to local absorption by the partial-covering absorber).

Download table as:  ASCIITypeset image

However, for both Models 1 and 2, the parameters that are allowed to vary between the high state and low state spectral models are still consistent within their 90% confidence intervals. Thus, although a flux-dependent spectral model is statistically motivated by the distribution of χ residual differences shown in Figure 12, the differences in spectral parameters between the high and low flux state are not substantial in magnitude and require better photon statistics to be properly constrained. These results are consistent with the low level of hardness/softness ratio variations as a function of count rate, and indicate that the fitting results of the flux-averaged spectrum reported in Section 4.1 are not substantially biased by the fact that the XMM-Newton and NuSTAR observations are not perfectly simultaneous.

As discussed in Section 4, the soft excess below 2 keV seen in the spectrum of IGR J18214-1318 can be accounted for either by introducing a blackbody component with properties typical of NS hot spots or a partial-covering absorber associated with the clumpy supergiant wind. In both models, the column density obscuring the whole binary system is measured to be ${N}_{{\rm{H}}}\approx 4\times {10}^{22}$ cm⁻², which can largely be ascribed to the interstellar medium and is consistent with the low column density measured by Rodriguez et al. (2009). The high column density measured by Tomsick et al. (2008), well in excess of the ISM value, is comparable to the ${N}_{{\rm{H}}}$ of the partial-covering absorber. Thus, the partial-covering model can naturally explain the observed variations in ${N}_{{\rm{H}}}$ as the result of changes in the density of clumps in the supergiant wind or how deeply embedded the compact object is in the stellar wind at different orbital phases.

However, these observed spectral variations are more difficult to explain using the blackbody model for the soft excess. Since spectra from the 2008 and 2009 soft X-ray observations of IGR J18214-1318 were fit with simple absorbed power-law models due to their low photon statistics, it is possible that variations in the strength of the blackbody emission could be incorrectly interpreted as ${N}_{{\rm{H}}}$ variations. We made fake 0.3–10 keV spectra with blackbody components of different strengths and fit them with simple absorbed power-law models to determine the effect that the blackbody emission alone can have on the measured ${N}_{{\rm{H}}}$, but found that it can only account for about 25% of the measured variations, which span the range of 3–12 × 10²² cm⁻². Although this does not rule out the presence of blackbody emission from this source, it suggests that there are periods of time when significant local obscuration is required to explain the observed ${N}_{{\rm{H}}}$ measurements. It is possible that both a partial-covering absorber and blackbody emission are present in this source but the photon statistics of these observations are insufficient to constrain the contribution of both components simultaneously.

In addition to the fact that the partial-covering absorber model provides a more natural explanation for the observed ${N}_{{\rm{H}}}$ variations, some unusual properties of the blackbody model make the partial-covering model the marginally preferred interpretation for the soft excess. Although the temperature (${kT}\sim 1.7$ keV) and size (${R}_{\mathrm{BB}}\approx 0.3$ km) of the blackbody emission region are typical for NS hot spots, emission from these hot spots is only seen in NS HMXBs with Be or main-sequence donors and ${L}_{{\rm{X}}}\lt {10}^{35}$ erg s⁻¹ (Reig et al. 2009; Bartlett et al. 2013). Furthermore, the flux of the blackbody component of IGR J18214-1318 is 73% ± 10% of the total 0.3–10 keV flux, which is much higher than the ∼30% hot spot blackbody flux fraction typically seen in other HMXBs (e.g., Mukherjee & Paul 2005; La Palombara & Mereghetti 2006, 2007; La Palombara et al. 2009). Some Sg HMXBs do exhibit a soft excess, but it is associated with blackbody emission with a lower temperature (${kT}\sim 0.2$ keV) from a larger area (${R}_{\mathrm{BB}}\sim 100$ km; Reig et al. 2009). Such emission is thought to arise from a cloud of diffuse photoionized plasma around the compact object associated with the supergiant wind; the photoionized plasma only absorbs photons at $\gtrsim 2\,\mathrm{keV}$ and produces significant emission at $\lesssim 1\,\mathrm{keV}$, resulting in a soft excess (Hickox et al. 2004; Szostek & Zdziarski 2008). However, the soft excess of IGR J18214-1318 has a higher blackbody temperature and smaller radius than is seen in Sg HMXBs. As a result, we believe that the partial-covering absorber is the most natural way of accounting for the soft excess, even if the blackbody model cannot be definitively dismissed.

Since neither pulsations nor cyclotron lines are detected in the XMM-Newton and NuSTAR data of IGR J18214-1318, we cannot definitively identify the compact object in this HMXB. Although we cannot rule out that this system hosts a BH, the exponential cutoff to its power-law spectrum with e-folding energy $\lt 25\,\mathrm{keV}$ argues in favor of an NS since BH HMXBs exhibit power-law spectra out to $\gtrsim 100\,\mathrm{keV}$ (Zdziarski 2000). Fitting the NuSTAR spectrum above 20 keV with a power-law model, we find that ${\rm{\Gamma }}=2.5\pm 0.2$. This soft photon index above the cutoff energy is typical for "normal" accreting pulsars, whereas anomalous X-ray pulsars (AXPs), which are thought to be magnetars, have Γ = 1–2 above 20 keV (Reig et al. 2012). The persistence of the hard X-ray emission from IGR J18214-1318 also disfavors a magnetar origin for this source, since most magnetars (all soft gamma-ray repeaters (SGRs) and many AXPs) exhibit bursting behavior (Olausen & Kaspi 2014 and references therein). Thus, the compact object in IGR J18214-1318 is most likely an accreting NS with a typical magnetic field strength (10${}^{12}\mbox{--}{10}^{13}$ G). The lack of detected pulsations with periods ≲1 hr can be explained by the geometry of the system (e.g., the nearly perfect alignment of the magnetic and rotational axes of the NS, the NS beam being narrow and not pointing toward Earth, or the NS beam being broad enough so as to wash out spin modulations) or a spin period >1–2 hr, longer than is typically seen in NS HMXBs.

IGR J18214-1318 is one of only a handful of HMXBs that may host an NS but for which sensitive periodicity searches have not discovered pulsations in the spin period range of 0.1–2000 s, which is typical for NSs in HMXBs. We refer to these HMXBs as non-pulsating, although it cannot be ruled out that they have pulsations with longer periods than are probed by current observations. Table 4 provides an overview of the properties of non-pulsating HMXBs, as well as HMXBs with long pulsation periods (${P}_{\mathrm{spin}}\gt 1\,\mathrm{hr}$) as a point of comparison.

Table 4.  Other Non-pulsating and Long Pulsation Period HMXBs

Source name	${P}_{\mathrm{spin}}$	${P}_{\mathrm{orbit}}$	Donor	Distance	LX (1035	${N}_{{\rm{H}}}$	Γ	${E}_{\mathrm{fold}}$	References
	(ks)	(days)	Type	(kpc)	erg cm−2 s−1)	(1022 cm−2)		(keV)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
Non-pulsating HMXBs
4U 1700-377	>50	3.41	O6.5 Iaf	1.9	0.4–40 (0.5–12 keV)	5–80	0.6–1.9	10–24	Gottwald et al. (1986), Heap & Corcoran (1992), Reynolds et al. (1999), Ankay et al. (2001), Boroson et al. (2003), van der Meer et al. (2005)
IGR J08262-3736	>5	⋯	OB V	6.1	0.2–0.4 (2–10 keV)	0.6–1.4	1.3–1.8	>70	Masetti et al. (2010), Bozzo et al. (2012)
IGR J16207-5129	>7	⋯	B1 Ia	6.1	0.2–5.0 (0.5–10 keV)	15–17	0.9–1.3	14–28	Nespoli et al. (2008), Tomsick et al. (2009), Ibarra et al. (2007), Bodaghee et al. (2010)
IGR J16318-4848	>10	⋯	sgB[e]	4.8	1.7 (20–50 keV) 66 (5–60 keV)	100–200	0.5–1.5	17–60	Walter et al. (2003), Filliatre & Chaty (2004), Barragán et al. (2009)
IGR J19140+0951	>2	13.55	B0.5 Ia	5.0	0.3–15 (2–20 keV)	<1–20	1.0–2.0	6–10	Hannikainen et al. (2007), Prat et al. (2008)
Long ${P}_{\mathrm{spin}}$ HMXBs
2S 0114+650	9.7	11.59	B1 Ia	4.5	1.5 (2–10 keV) 4.3 (3–20 keV)	2–5	0.8–1.2	13–18	Reig et al. (1996), Corbet et al. (1999), Farrell et al. (2008)
4U 2206+54	5.56	9.56	O9.5 V	2.6	0.35–3.5 (1–12 keV)	0.9–4.7	0.9–1.8	10–29	Blay et al. (2005), Ribó et al. (2006), Reig et al. (2009)

Notes. (2) Pulsation period. For non-pulsating HMXBs, the value provided is the longest period that has been ruled out by periodicity searches. It is possible that these sources have pulsation periods longer than these values. (3) Orbital period of the binary (for sources for which it has been measured). (4) Spectral type of the donor star. (5) Distance estimate. Typical uncertainty is ±2 kpc. (6) Range of observed X-ray luminosities for the source. The energy band is written in parentheses. (7)–(9) Range of observed spectral parameters based on different observations of the source.

Download table as:  ASCIITypeset image

As can be seen in Table 4, all but one of the non-pulsating HMXBs exhibit exponential cutoffs to their power-law spectra with e-folding energies $\lesssim 30\,\mathrm{keV}$, suggesting that they likely harbor NSs. The one exception is IGR J08262-3736, for which no exponential cutoff is measured below 70 keV (Bozzo et al. 2012). The lack of an exponential cutoff in the spectrum of IGR J08262-3736 makes it the most plausible BH candidate of all the non-pulsating HMXBs. This HMXB also differs from all the other non-pulsating HMXBs, including IGR J18214-1318, in that its donor star is a main-sequence OB star (Masetti et al. 2010).

Of the four non-pulsating Sg HMXBs, IGR J16207-5129 most resembles IGR J18214-1318. In addition to having a similar X-ray luminosity, photon index, and exponential cutoff energy, this source also exhibits a soft excess, which Tomsick et al. (2009) argue most likely results from partial obscuration by the stellar wind with ${N}_{{\rm{H}}}\approx {10}^{23}$ cm⁻². The PDS of IGR J16207-5129, like the PDS of IGR J18214-1318, shows significant red noise below 0.01 Hz, which is well-fit by a simple power law; the slope of the red noise in IGR J16207-5129 is steeper, having $\alpha =1.76\pm 0.05$ (Tomsick et al. 2009), but still falling within the range that is typical for frequencies above the pulsation frequency, suggesting, but not proving, that the NS in this HMXB may have a period longer than the ∼2 hr limit probed by the data.

Another non-pulsating Sg HMXB that is similar to IGR J18214-1318 is IGR J19140+0951, except for the fact that IGR J19140+0951 exhibits a soft excess modeled by a lower temperature (kT = 0.3 keV) blackbody, which is attributed to a shock formed between the ionized gas around the NS and the stellar wind (Prat et al. 2008). The ${N}_{{\rm{H}}}$ variations measured in the spectrum of IGR J19140+0951 are related to its orbital phase and are even larger than those observed for IGR J18214-1318.

Apart from having similar photon indices and exponential cutoff energies indicative of NS HMXBs, the other two non-pulsating Sg HMXBs have more extreme properties than IGR J18214-1318. IGR J16318-4848 is the HMXB with the highest measured local column density (${N}_{{\rm{H}}}\sim {10}^{24}$ cm⁻²; Walter et al. 2003), which is due to its sgB[e] donor star (Filliatre & Chaty 2004). The relatively slow 400 km s⁻¹ wind of this donor star helps maintain a high mass accretion rate onto the compact object, resulting in a slightly higher X-ray luminosity (${L}_{{\rm{X}}}\sim {10}^{36}$ erg s⁻¹) than in other Sg HMXBs (Barragán et al. 2009). 4U 1700-377 displays the highest levels of variability of any of the non-pulsating Sg HMXBs; its flux, even at hard X-ray energies, varies by factors >100, and the pattern of spectral variations as a function of luminosity provides further evidence that 4U 1700-377 likely hosts an NS (Seifina et al. 2016). However, the mass of the compact object in this system has been measured to be 2.44 ± 0.27 M_⊙ (Clark et al. 2002), placing it among the most massive NSs observed (Ozel & Freire 2016). For 4U 1700-377 and IGR J16318-4848, the existence of pulsations is constrained over a wider range of periods than for other sources, excluding all periods <50 ks and <10 ks, respectively.

It is possible that some of the non-pulsating HMXBs, especially those for which periodicity searches do not extend beyond 1–2 hr, possess slowly spinning NSs. Currently, two HMXBs, 4U 2206+54 and 2S 0114+650, are known to have pulsation periods longer than an hour (Corbet et al. 1999; Reig et al. 2009). The photon indices and cutoff energies of both sources are similar to IGR J18214-1318 (Farrell et al. 2008; Reig et al. 2009), although 2S 0114+650 is a closer analogue since it has a supergiant donor (Reig et al. 1996). Li & van den Heuvel (1999) proposed that the slow spin of the NS in 2S 0114+650 indicates that it was born as a magnetar with B ≳ 10¹⁴ G, was slowed down efficiently by the propeller effect before its magnetic field significantly decayed to its current expected value of ∼10¹² G. In the case of 4U 2206+54, magnetorotational models, which can account for the NS's spin and spin-down rate, require magnetic fields strengths between 5 × 10¹³ and 3 × 10¹⁵ G (Ikhsanov & Beskrovnaya 2010); thus, it is possible that 4U 2206+54 currently contains a magnetar and would evolve into a system like 2S 0114+650. Neither source exhibits clear cyclotron absorption features, although low-significance detections of such features have been claimed by some authors, which would indicate that both HMXBs host NSs with typical strength magnetic fields (Torrejón et al. 2004; Bonning & Falanga 2005).

NSs with long spin periods have also been discovered in some symbiotic X-ray binaries (SyXBs). The two SyXBs with the longest spin periods, IGR J16358-4726 and 4U 1954+319, have spin periods of ∼1.6 hr and ∼5.3 hr, respectively (Kouveliotou et al. 2003; Patel et al. 2004; Corbet et al. 2006; Marcu et al. 2011). Both of these sources have spectra that are very similar to those of long spin-period HMXBs and IGR J18214-1318 (Patel et al. 2007; Enoto et al. 2014), but they have giant M-type stellar companions rather than high-mass donors. Although it has been suggested that the long spin-period SyXBs may host magnetars or magnetar descendants just like long spin-period HMXBs (Patel et al. 2007; Enoto et al. 2014), models that assume quasi-spherical wind accretion for SyXBs rather than disk accretion do not require magnetar-strength fields to explain their timing properties (Lü et al. 2012).

In summary, in addition to IGR J18214-1318, there are five non-pulsating HMXBs, four of which have exponential cutoff energies <30 keV, suggesting they most likely harbor NSs. Like IGR J18214-1318, all four of these HMXBs have supergiant donor stars and resemble the Sg HMXB 2S 0114+650, which hosts an NS with a 2.7 hr pulsation period thought to have been born as a magnetar. Population synthesis models predict that 8%–9% of all NSs are born as magnetars, and that only ∼2% of NSs in binaries are magnetars; these models predict that an even a smaller percentage of magnetars would be part of an X-ray binary because many of them are produced from the secondary rather than the primary (Popov & Prokhorov 2006). About 100 HMXBs have been discovered in the Galaxy, and only about 60 of them are known to host pulsars (Bird et al. 2016), so with the discovery of 2S 0114+650 and 4U 2206+54, both of which may host magnetars (or former magnetars), the observed number of magnetar HMXBs already agrees with theoretical expectations. Thus, if future observations reveal that several of the "non-pulsating" HMXBs actually host long-period pulsars, it could imply that either current models of magnetar origins or models of the spin evolution of NSs in binaries need to be revised.