An XMM-Newton and NuSTAR study of IGR J18214-1318: a non-pulsating high-mass X-ray binary with a neutron star

TL;DR

This study uses XMM-Newton and NuSTAR observations to analyze IGR J18214-1318, a high-mass X-ray binary, suggesting its compact object is likely a neutron star despite no detected pulsations.

Contribution

It provides the first detailed spectral analysis of IGR J18214-1318, proposing the neutron star nature based on spectral features and variability, despite the absence of pulsation detection.

Findings

Spectral cutoff indicates a neutron star primary

Presence of Fe Kα line and soft excess in spectrum

No pulsations detected within the observational limits

Abstract

IGR J18214-1318, a Galactic source discovered by the International Gamma-Ray Astrophysics Laboratory, is a high-mass X-ray binary (HMXB) with a supergiant O-type stellar donor. We report on the XMM-Newton and NuSTAR observations that were undertaken to determine the nature of the compact object in this system. This source exhibits high levels of aperiodic variability, but no periodic pulsations are detected with a 90% confidence upper limit of 2% fractional rms between 0.00003-88 Hz, a frequency range that includes the typical pulse periods of neutron stars (NSs) in HMXBs (0.1-10$^3$ s). Although the lack of pulsations prevents us from definitively identifying the compact object in IGR J18214-1318, the presence of an exponential cutoff with e-folding energy $\lesssim30$ keV in its 0.3-79 keV spectrum strongly suggests that the compact object is an NS. The X-ray spectrum also shows a Fe K$\alpha$ emission line and a soft excess, which can be accounted for by either a partial-covering absorber with $N_{\mathrm{H}}\approx10^{23}$ cm$^{-2}$ which could be due to the inhomogeneous supergiant wind, or a blackbody component with $kT=1.74^{+0.04}_{-0.05}$ keV and $R_{BB}\approx0.3$ km, which may originate from NS hot spots. Although neither explanation for the soft excess can be excluded, the former is more consistent with the properties observed in other supergiant HMXBs. We compare IGR J18214-1318 to other HMXBs that lack pulsations or have long pulsation periods beyond the range covered by our observations.

Tables

Table 1. Table 1 Observations of IGR J18214-1318

Telescope	Observation ID	Start Time	Duration	Instrument/	Exposure
		(UTC)	(ks)	Detector	(ks)
(1)	(2)	(3)	(4)	(5)	(6)
NuSTAR	3000114002	2014 Sep 18 01:16:07	49.8	FPMA	26.0
				FPMB	25.8
XMM-Newton	0741470201	2014 Sep 18 02:34:26	26.9	EPIC pn	18.5
				MOS1	25.9
				MOS2	25.8

Table 2. Table 2 Best-fit Spectral Parameters

	Model 1	Model 2
	tbabs*(bbody+	tbabs*pcfabs*
	powerlaw*highecut	(powerlaw*highecut
	+gaussian)	+gaussian)
$N_H$ (${10}^{22}$ cm-2)	4.2${}_{-0.2}{}^{}{}_{}{}^{+0.3}$	4.3$\pm$0.6
$k{T}_{\mathrm{BB}}$ / $N_{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}$
$\Gamma$	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$\pm$0.02
${\sigma }_{\mathrm{line}}$ (eV)
EWline (eV)	53${}_{-21}{}^{}{}_{}{}^{+16}$	57${}_{-14}{}^{}{}_{}{}^{+16}$
$C_{\mathrm{MOS1},2}$	1.04$\pm$0.02	1.04$\pm$0.02
$C_{\mathrm{FPMA}}$	1.25$\pm$0.02	1.24$\pm$0.02
$C_{\mathrm{FPMB}}$	1.32$\pm$0.03	1.31$\pm$0.03
${\chi }_{\nu }^2$/dof	1.10/578	1.10/578

Table 3. 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_H$ (${10}^{22}$ cm-2)		4.4${}_{-0.4}{}^{}{}_{}{}^{+0.6}$	3.6${}_{-0.4}{}^{}{}_{}{}^{+0.7}$	4.0$\pm$1.1
$k{T}_{\mathrm{BB}}$ (keV) / $N_{H,\mathrm{partial}}$ (${10}^{22}$ cm-2)	1.74${}_{-0.18}{}^{}{}_{}{}^{+0.06}$				7$\pm 2$	11${}_{-3}{}^{}{}_{}{}^{+4}$
BB norm. (10-11 erg cm-2 s-1) / Cov. frac.		1.07$\pm 0.14$	0.85$\pm 0.15$		0.9$\pm 0.1$	0.8$\pm 0.1$
$\Gamma$	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$\pm 0.04$	2.80${}_{-0.21}{}^{}{}_{}{}^{+0.24}$		1.10$\pm 0.04$	2.63${}_{-0.11}{}^{}{}_{}{}^{+0.12}$
$C_{\mathrm{FPMB}}$		1.10$\pm 0.05$	2.99${}_{-0.23}{}^{}{}_{}{}^{+0.25}$		1.14$\pm 0.05$	2.81${}_{-0.12}{}^{}{}_{}{}^{+0.13}$
${\chi }_{\nu }^2$/dof	1.01/481			1.02/481

Table 4. Table 4 Other non-pulsating and long pulsation period HMXBs

Source name	$P_{\mathrm{spin}}$	$P_{\mathrm{orbit}}$	Donor	Distance	$L_X$ (1035	$N_H$	$\Gamma$	$E_{\mathrm{fold}}$	References
	(ks)	(days)	type	(kpc)	erg cm-2 s-1)	(${10}^{22}$ 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)

Equations

$$\mathrm{rms}=\sqrt{\frac{P_{\mathrm{UL}}-2}{CR}{\left(\frac{S+B}{S}\right)}^2\Delta \nu }$$

Full text

††slugcomment: Accepted to ApJ: April 27, 2017

I. Introduction

High-mass X-ray binaries (HMXBs) inform our understanding of the evolution of massive stars, which is still subject to significant uncertainties (Smith, 2014). Studying the accreting neutron stars (NSs) and black holes (BHs) in these systems offers a special tool to probe the strength and clumping of the stellar winds of their massive companions. Moreover, comparing the properties of HMXB populations to predictions of population synthesis models helps constrain theoretical models of stellar mass loss, mass transfer episodes in massive binaries, and the natal kicks received by compact objects during supernova explosions (e.g., Negueruela et al. 2008; Linden et al. 2010; Fragos et al. 2013; Grudzinska et al. 2015). Constraining the ratio of NS to BH HMXBs and whether this ratio varies with binary properties (e.g., donor spectral type and metallicity) can shed light on the net mass loss experienced by a high-mass star due to its stellar wind, binary interactions, and supernova explosion (Dray 2006; Muno 2007; Belczynski & Ziolkowski 2009) as well as improve our estimates of the relative fractions of different double compact binaries expected to descend from HMXBs (Postnov & Yungelson, 2014).

HMXBs hosting BHs exhibit different spectral properties from those hosting NSs. NS HMXBs typically have hard power-law spectra with exponential cutoffs, with e-folding energies typically $\lesssim 30$ keV (Coburn et al., 2002), whereas BH HMXBs exhibit power-law cutoffs around 50–100 keV in their hard states and $\Gamma\sim 2$ power-law tails extending to MeV energies in their soft states (Grove et al. 1998; Zdziarski 2000). Thus, the presence of an exponential cutoff below 20 keV in the X-ray spectrum of an HMXB is a strong indication that it harbors an NS; however, only the detection of X-ray pulsations or cyclotron line features constitute definitive proof of the presence of an NS Most X-ray pulsars in HMXBs have spin periods between $\sim 0.1$ and $\sim 10^{3}$ s (Corbet 1986; Chaty 2013), although a couple of longer-period pulsars have been discovered (Reig et al. 2009; Corbet et al. 1999). In addition to NS X-ray pulsations, the X-ray light curves of some HMXBs can exhibit orbital or superorbital modulations with typical periods of a few hours to a few hundred days (Corbet et al. 2006; Corbet & Krimm 2013). The cyclotron lines that have been observed in some HMXB spectra have energies between 10 and 80 keV, corresponding to magnetic field strengths of a few $10^{12}$ G (Coburn et al. 2002; Pottschmidt et al. 2005; Caballero et al. 2007; Doroshenko et al. 2010; Caballero & Wilms 2012; Tsygankov et al. 2012; Fürst et al. 2014; Tendulkar et al. 2014; Yamamoto et al. 2014; Bellm et al. 2014). An HMXB with an unbroken power-law spectrum extending beyond 50 keV that does not show X-ray pulsations can be considered a BH candidate, but confirming the BH nature of the compact object requires a dynamical mass measurement in excess of 2–3 $M_{\odot}$, the maximum theoretically expected NS mass (Lattimer, 2012); however, such measurements can be challenging to obtain.

Since its launch in 2002, the International Gamma-Ray Astrophysics Laboratory (INTEGRAL) has discovered a large number of sources that were given ”IGR” source names, including 40 new HMXBs in the 4th IBIS/ISGRI catalog (Bird et al., 2016). X-ray pulsations have only been detected from about a quarter of the IGR HMXBs, and the nature of the compact object in the remaining systems is undetermined since many of them are relatively recent discoveries and have not yet been well-studied. Identifying the nature of compact objects in IGR HMXBs is of special interest since they are demographically different from the population of HMXBs discovered prior to INTEGRAL, consisting of roughly equal numbers of HMXBs with Be and supergiant (Sg) stellar companions, whereas Be HMXBs dominate the pre-INTEGRAL population (Walter et al., 2015). Many of the IGR Sg HMXBs have unusual properties, exhibiting high levels of obscuration ($N_{\mathrm{H}}\sim 10^{23}-10^{24}$ cm*-2*) local to the source (Walter et al. 2006; Chaty et al. 2008) or extreme flaring behavior characterized by hard X-ray flux variations of several orders of magnitude on timescales of a few hours (Negueruela et al. 2006; Sguera et al. 2006). The discovery of these Sg HMXBs was made possible by the greater sensitivity and higher cadence of the INTEGRAL Galactic Plane survey compared to previous hard X-ray missions.

A few of the flaring IGR Sg HMXBs, known as supergiant fast X-ray transients (SFXTs), are known to host NSs based on the detection of X-ray pulsations (Romano et al., 2014) or cyclotron lines (Bhalerao et al., 2015), and several of the models proposed to explain SFXT behavior depend on the presence of an NS magnetosphere (the propeller effect, Grebenev & Sunyaev 2007; magnetic gating, Bozzo et al. 2008; quasi-spherical settling accretion of hot plasma shells, Shakura et al. 2014), suggesting that all SFXTs may host NSs. However, the nature of compact objects in many IGR HMXBs that are not SFXTs remains unknown, and as several Sg HMXBs are known to harbor BHs (e.g., Cyg X-1, M33 X-7, LMC X-1, LMC X-3) compared to only one known Be-BH binary (Casares et al., 2014), IGR Sg HMXBs constitute a particularly promising group to search for BHs.

IGR J18214-1318 is one of the Sg HMXBs lacking a clear compact object identification. This source was first reported in the second IBIS/ISGRI catalog (Bird et al., 2006) and detected consistently by INTEGRAL with a flux of 1–2 mCrab in the 20-40 keV band (Krivonos et al. 2012; Bird et al. 2016). The source was localized with arcsecond precision to R.A. = $18^{\mathrm{h}}21^{\mathrm{m}}19.76^{\mathrm{s}}$, decl. = $-13^{\circ}18^{\prime}38.9^{\prime\prime}$ through a Chandra observation (Tomsick et al., 2008). The localization of this source permitted its association with an optical counterpart, which is a high-mass star of most likely spectral type O9 I (Butler et al., 2009), thus securing the identification of IGR J18214-1318 as a Sg HMXB. Its Chandra spectrum is well-fit by an absorbed power law with $N_{\mathrm{H}}=1.2\pm 0.3\times 10^{23}$ cm*-2* and $\Gamma=0.7^{+0.6}_{-0.5}$. Later Swift observations of this source measured a similar photon index ($\Gamma=0.4\pm 0.2$) but a much lower absorbing column density of $N_{\mathrm{H}}=3.5^{+0.8}_{-0.5}\times 10^{22}$ cm*-2*, which is consistent with the Galactic $N_{\mathrm{H}}$ integrated along the line of sight (3.1$\times 10^{22}$ cm*-2*)111The Galactic $N_{\mathrm{H}}$ along the line of sight to IGR J18214-1318 is calculated as the sum of the HI contribution measured from the Leiden/Argentine/Bonn (LAB) survey of HI (Kalberla et al., 2005) and the H2 contribution estimated from the MWA CO survey (Bronfman et al., 1989).. Although the Chandra data suffered from photon pile-up, since the photon indices derived by Chandra and Swift are so similar, it is unlikely that the large difference in the derived $N_{\mathrm{H}}$ values is simply a result of the photon pile-up in Chandra. Thus, the large variability of $N_{\mathrm{H}}$ is likely real and associated with material local to the source; similar $N_{\mathrm{H}}$ variations have been seen in other Sg HMXBs (e.g, IGR J19140+0951; Prat et al. 2008). The hard power-law index measured in the soft X-ray band suggests that the compact object in IGR J18214-1318 is more likely to be an NS than a BH, but it does not constitute strong or definitive evidence.

Therefore, we observed IGR J18214-1318 with the Nuclear Spectroscopic Telescope Array (NuSTAR) and XMM-Newton (§II) to better constrain the nature of the compact object in this HMXB. NuSTAR and XMM-Newton are ideally suited for this study, because their instruments have the fast temporal resolution required to search for X-ray pulsations and their combined broadband X-ray spectral coverage from 0.3-79 keV permits the measurement of cyclotron lines and hard X-ray cutoffs, which may be present in the HMXB power-law spectrum. Spectral analysis of this data, resulting in the detection of a high-energy cutoff, is described in §IV, and timing analysis ruling out the existence of a pulse period $\lesssim 1$ hr is presented in §III. In §V, we discuss IGR J18214-1318 in the context of other non-pulsating and long pulse period HMXBs.

II. Observations and Data Reduction

NuSTAR and XMM-Newton observed IGR J18214-1318 on 2014 September 18. Observation details are provided in Table 1. The duration of the NuSTAR observation is about twice as long as that of XMM-Newton because Earth occultations reduce the effective exposure of NuSTAR. The exposure times of the two NuSTAR focal plane modules (FPM) and the two XMM Metal Oxide Semi-conductor (MOS) CCD arrays of the European Photon Imaging Camera (EPIC) are roughly equal (26 ks), while the exposure time of the EPIC pn CCD is significantly lower (18.5 ks) due to its higher dead-time fraction.

II.1. NuSTAR

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*′′* 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*′′* and an outer radius of 250*′′. Within 60′′* 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*′′* radius (corresponding to the PSF encircled energy fraction of 75%) to limit background contamination while still ensuring good photon statistics ($\sim$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*-1* for FPMA and FPMB combined.

II.2. XMM-Newton

We made use of data from the EPIC pn and MOS instruments on XMM-Newton, which provide 15*′′* 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 similar energy resolution in the 0.3–10 keV band but have 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 (SAS) 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, 222See http://xmm.esac.esa.int/sas/current/documentation/threads/, 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^{\prime\prime}$ 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^{\prime}$ 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*-1* with EPIC pn and 0.8 counts s*-1* for the two MOS cameras combined. We checked the observations for pile-up and found that it was not an issue.

III. Timing Analysis

IGR J18214-1318 exhibits a high level of variability in its X-ray light curve, as can be seen in Figure 2. This strong aperiodic variability is common in HMXBs (Belloni & Hasinger, 1990), and is attributed to variations in the accretion rate resulting from density perturbations due to either disk instabilities or magnetohydrodynamic turbulence (e.g. Shakura & Sunyaev 1976; Hoshino & Takeshima 1993; Revnivtsev et al. 2009). In order to determine whether this HMXB hosts an NS, we searched for periodic pulsations in this noisy light curve through analysis of its power spectrum. We performed this pulsation search using both NuSTAR and XMM-Newton EPIC pn data, since the EPIC pn camera has the highest temporal resolution and effective area of the XMM instruments. Although NuSTAR’s temporal resolution is better than XMM’s, the Earth occultations that NuSTAR experiences during its orbit create large gaps in its light curves, which in turn introduce additional noise in the power spectrum at low frequencies. Furthermore, since pulsars in Sg HMXBs tend to have periods $\gtrsim$1 s (Corbet 1986; Skinner et al. 1982), the 5.7 ms resolution of XMM-Newton should be sufficient for detecting possible pulsations in IGR J18214-1318. We first describe our analysis of the XMM EPIC pn power spectrum, and then compare the results to those obtained using NuSTAR data.

We corrected the arrival time of each event detected by the XMM-Newton EPIC pn instrument to the solar system barycenter, and used this corrected event list to make a light curve in the 0.3–12 keV band with the maximum possible time resolution of 5.6718 ms. We then used the XRONOS tool powspec to produce a Leahy-normalized power spectrum (Leahy et al., 1983) of this light curve, shown in Figure 3. The power spectrum spans frequencies from 3.7$\times 10^{-5}$ Hz (based on the 27 ks duration of the observation) to 88.1 Hz (the Nyquist frequency). In a Leahy-normalized power spectrum, Poissonian noise results in power being distributed as a $\chi^{2}$ probability distribution with two degrees of freedom (dof); we used this distribution and the number of trials (which is equal to the number of frequency bins) to calculate the 90% confidence detection threshold for this power spectrum as 35.1, shown by the dashed lines in Figure 3. This detection threshold is only exceeded at frequencies $<5.06\times 10^{-3}$ Hz, but significant red noise at low frequencies suggests we should be cautious in ascribing this excess power to periodic pulsations. Thus, in order to account for the red noise present at frequencies $<$0.1 Hz, we analyzed the power spectrum above and below 0.1 Hz separately.

Above 0.1 Hz, the maximum Leahy power ($P_{\mathrm{max}}$) measured is 30.9, which is below the 90% confidence threshold but can be used to calculate an upper limit on the strength of a periodic signal. As derived by van der Klis (1989), the 90% confidence upper limit of the Leahy power ($P_{\mathrm{UL}}$) is given by $P_{\mathrm{UL}}~{}=~{}P_{\mathrm{max}}~{}-~{}P_{\mathrm{exceed}}$, where $P_{\mathrm{exceed}}$ is the power level exceeded by 90% of the frequency bins. In our case, $P_{\mathrm{exceed}}$ is 0.2, which implies that $P_{\mathrm{UL}}$ is 30.7. This upper limit on the Leahy power can be converted into an upper limit on the source fractional rms variability for a periodic signal using the following formula:

[TABLE]

where $CR$ is the mean (source plus background) count rate, $S$ represents the net source counts, $B$ represents the estimated background counts in the source region, and $\Delta\nu$ is the width of the frequency bin. Thus, the 90% upper limit on the rms noise level for a periodic signal between 0.1 and 88 Hz is $<2.2$%.

In order to search for periodic signals at low frequencies, we first characterized the red noise below 0.1 Hz. We produced a 0.3–12 keV light curve with 5 s resolution, and used it to make a rms-normalized power density spectrum (PDS) with a Nyquist frequency of 0.1 Hz and a minimum frequency of 9.8$\times 10^{-5}$ Hz, since we averaged together the PDS made from three time intervals of $\sim 10$ ks. The rebinned PDS from which the Poisson noise level of 2.0 has been subtracted is shown in Figure 4. The PDS is dominated by red noise below 0.02 Hz, so we fit the PDS below this frequency with a power-law model ($P=A(\nu/\mathrm{1Hz})^{-\alpha}$) using Whittle statistics (Whittle 1953; Whittle 1957; Vaughan 2010). The best-fit parameters are $\alpha=1.40^{+0.02}_{-0.20}$ and $A=0.0018^{+0.006}_{-0.0001}$, where the quoted errors correspond to the 90% confidence intervals. The integrated source fractional rms for frequencies between 10*-3* to 0.1 Hz is $24$%$\pm 1$%. The rms PDS slope and integrated fractional rms measured using 3–12 keV XMM light curves are consistent with the values measured from the 0.3–12 keV data to better than 1$\sigma$ confidence. Fitting a power-law model plus a constant to the Leahy-normalized power spectrum from 9.8$\times 10^{-5}$ to 88 Hz results also results in a consistent power-law slope. We also tried fitting a broken power-law model to the rms-normalized and the Leahy-normalized power spectra, but in both cases, the best-fitting break frequency was poorly constrained and exceeded the maximum frequency of the power spectrum. Thus, we do not find a significant break in the power spectrum between 0.0003 and 88 Hz.

In order to properly subtract the red noise component from the Leahy-normalized power spectrum in the 0.00003–0.1 Hz range, we follow the procedure described in van der Klis (1989). The Leahy power spectrum was multiplied by a factor of 2 and then divided by the (appropriately renormalized) best-fit power-law model; we checked that the resulting power spectrum, shown in Figure 5, follows a $\chi^{2}$ distribution with 2 dof, a requirement for applying the method described in van der Klis (1989) to power spectra with red noise. Taking into account the number of trials, we calculated the 90% confidence detection threshold to be 21.2. As can be seen in Figure 5, no frequency bin exceeds this power level. Given that $P_{\mathrm{max}}$ in this low-frequency range is 17.0 and $P_{\mathrm{exceed}}$ is 0.2, $P_{\mathrm{UL}}$ is 16.8, which corresponds to an upper limit on the source fractional rms of 1.7%. We verified that varying the red noise slope and normalization within their 90% confidence intervals does not significantly affect these results.

In some binaries, a periodic signal may be difficult to detect because it may be spread out in frequency space due to orbital motion. However, HMXBs tend to have orbital periods of a few to hundreds of days, and thus our observations of IGR J18214-1318 are short enough that orbital modulations of the pulsation frequency should not be significant. Thus, given the stringent upper limits on the fractional rms for periodic signals in the 0.00003–88 Hz frequency range, pulsations with periods $\lesssim$1 hr are strongly ruled out.

We performed the same timing analysis with NuSTAR data in the 3–12 and 12–30 keV bands. All NuSTAR photon arrival times were converted to barycentric dynamical time (TDB). After being corrected for thermal drift of the on-board clock, the NuSTAR time resolution is $\sim 2$ ms rms, and its absolute accuracy is known to be better than 3 ms (Mori et al., 2014). Thus, we used the NuSTAR data to search for spin periods as short as $\sim$1 ms by binning the light curves by 1/2048 s (488 $\mu$s). We produced Leahy-normalized power spectra from FPMA and FPMB light curves using both the XRONOS tool powspec and the power spectrum tools developed by M. Bachetti, 333Tools can be found at https://bitbucket.org/mbachett/maltpynt and are described in Bachetti et al. (2015). but we did not find any significant peaks in the power spectrum below 1024 Hz. Above 0.1 Hz, the 90% confidence upper limit on the source fractional rms is 3.4% in the 3–12 keV band and 6.4% in the 12–30 keV band.

We produced an rms-normalized, noise-subtracted PDS and fit the red-noise-dominated continuum in the $6\times 10^{-5}$ to 0.02 Hz frequency range with a power-law model. Jointly fitting the PDS produced from FPMA and FPMB data yields $\alpha=1.32^{+0.06}_{-0.10}$ in the 3–12 keV band and $\alpha=1.46^{+0.07}_{-0.12}$ in the 12–30 keV band. These measured power-law slopes are consistent with the values derived from XMM-Newton. Using the PDS power-law fits to “normalize” the red noise continuum in the Leahy power spectrum below 0.1 Hz, we checked that the resulting power spectrum follows a $\chi^{2}$ distribution with 2 dof, and then calculated that the 90% confidence upper limit on the source fractional rms is 2.7% in the 3–12 keV band and 6.0% in the 12–30 keV band. The integrated source fractional rms between 10*-3* to 0.1 Hz is 31%$\pm$2% and 23%$\pm$8% in the 3–12 and 12–30 keV bands, respectively. Thus, the NuSTAR 3–12 keV integrated fractional rms is higher at 3$\sigma$ confidence than that measured in the XMM-Newton 3–12 keV band; this difference may be partly attributed to the XMM-Newton and NuSTAR observations not being fully coincident in time and the additional artificial noise injected into the NuSTAR power spectrum by the light curve gaps due to Earth occultations. Due to the large errors bars of the integrated fractional rms measured in the NuSTAR 12–30 keV band, this value is consistent at 1$\sigma$ confidence with both the integrated rms measured by XMM-Newton and the NuSTAR 3–12 keV band.

Overall, the results from our XMM-Newton and NuSTAR timing analyses are in agreement, ruling out the presence of pulsations with periods shorter than about an hour. The integrated source fractional rms values measured between $10^{-3}$ and 0.1 Hz by both telescopes are within the typical range of 10-30% seen in HMXBs (Belloni & Hasinger, 1990). For most accreting X-ray pulsars, the red noise power-law index is $\alpha=1.4-2.0$ at frequencies higher than the pulsation frequency and $\alpha=0-1.0$ at lower frequencies (Hoshino & Takeshima, 1993). Thus, the red noise power-law slope of IGR J18214-1318 is similar to the slopes observed above the pulsation frequency in X-ray pulsars. This fact, combined with the lack of a frequency break in the red noise continuum, suggests, but does not prove, that the pulsation frequency in IGR J18214-1318 may be lower than the range probed by our data.

IV. Spectral Analysis

Since IGR J18214-1318 exhibits strong variability, when performing spectral analysis it is important to consider whether the source spectrum varies with source brightness. First, we checked for substantial spectral variations correlated with source flux by producing light curves in soft, medium, and hard energy bands (0.3–3, 3–12, and 12–30 keV, respectively), and using them to calculate hardness/softness ratios as a function of time, which are shown in Figure 6. With the NuSTAR data, we calculate the hardness ratio as the count rate in the hard energy band divided by the rate in the medium band. With the XMM-Newton data, we calculate the softness ratio as the count rate in the soft energy band divided by the rate in the medium energy band. The fact that the XMM-Newton softness ratio and the NuSTAR hardness ratio have the same energy band in the denominator helps to visualize how the soft and hard energy ends of the X-ray spectrum vary relative to one another. The average ratios in both the XMM-Newton and NuSTAR bands do not vary significantly on timescales of $3-4$ ks, remaining consistent within 1$\sigma$ confidence intervals. However, there are indications of ratio variability on shorter timescales, as several XMM-Newton and NuSTAR ratios (5 of 265 and 10 of 287, respectively) measured in 100s intervals differ from the mean ratio by $>2.7\sigma$, even though only one ratio measurement per instrument is statistically expected to differ by this amount.

In order to establish whether there is a significant correlation between the source brightness and X-ray ratios, we grouped together the individual measurements from 100 s light-curve fragments with similar count rates, creating a total of five count rate bins and evaluating the mean X-ray hardness/softness ratio in each bin. The standard deviation of the distribution of ratio values within each count rate bin differs by $<15$% from the average 1$\sigma$ uncertainty of the ratio measurements in that bin, indicating that the apparent scatter in ratios within each count rate bin is primarily statistical, not intrinsic, in nature; thus, averaging the ratio values within a given count rate bin reduces the statistical error, which for each individual ratio measurement is $\pm 0.1$ on average. We calculate the error on the mean as the quadrature sum of the individual measurement uncertainties divided by the total number of measurements averaged together.

As shown in Figure 7, in both the XMM-Newton and NuSTAR bands, the mean X-ray ratio increases slightly with increasing count rate, as is observed in some other Sg HMXBs (e.g., IGR J16207$-$5129 and 4U 2206+54; Tomsick et al. 2009; Wang 2013). Fitting a constant value to the mean ratios in Figure 7 using $\chi^{2}$ minimization, we find that the null hypothesis that the ratios do not vary with count rate is rejected at 98.4% confidence in the XMM-Newton band and at 99.9% confidence in the NuSTAR band. While the variations of X-ray ratios with count rate are statistically significant, the spectral variations of IGR J18214-1318 are likely to be small given that the mean X-ray ratios increase by only 25% as the count rate increases by almost an order of magnitude; therefore, it is unlikely that significant biases would be introduced by fitting the average source spectrum derived from all available data. Nonetheless, verifying whether the flux-dependent spectral variations are indeed small and determining the cause of variations, if they exist, can be valuable for understanding the origin of variability in IGR J18214-1318. Therefore, we made one spectrum, which we refer to as the flux-averaged spectrum, combining all XMM-Newton and NuSTAR data to achieve maximum sensitivity to detect potentially weak features such as emission lines or cyclotron absorption lines, and we investigated the origin and significance of spectral variations by making a low flux and a high flux spectrum based on simultaneous intervals of XMM-Newton and NuSTAR observations.

IV.1. The flux-averaged spectrum

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 §II. The spectra were rebinned with the requirement that the signal significance in each bin be $\geq 10$, except for the highest energy bin which was required to have a significance $\geq 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$\sigma$ 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^{2}_{\nu}$ = 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.

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$\alpha$ line emission at $6.40\pm 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.

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 the 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.21}_{-0.18}$ and line energy $E_{\mathrm{cyc}}=27^{+4}_{-6}$ keV; the inclusion of cyclabs only reduced the chi-squared value of the fit by 5 and left $\chi^{2}_{\nu}$ 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^{+4}_{-7}$ 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^{2}_{\nu}$ 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 datasets, 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 dataset 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 ($\Delta\chi^{2}$) between the two model fits was calculated. The maximum value of $\Delta\chi^{2}$ from these simulations was 19.3, slightly higher than the observed value. Based on the distribution of $\Delta\chi^{2}$ from our simulations, we estimate there is roughly a 0.001% chance of measuring the observed value of $\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 $<25$ 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$ 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\times 10^{22}$ cm*-2*, 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 $\sim 10^{23}$ cm*-2*. 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.02}_{-0.05}\times 10^{-11}$ erg cm*-2* s*-1* in the XMM-Newton EPIC pn observations and 2.18${}^{+0.03}_{-0.10}\times 10^{-11}$ erg cm*-2* s*-1* in the NuSTAR observations (FPMA and FPMB averaged). At a distance of 9–10 kpc, these fluxes correspond to unabsorbed luminosities of $\approx 1-2\times 10^{35}$ erg s*-1*. 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.

IV.2. Spectral variations with brightness

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 100s 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*-1*. 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 $\geq 7$, except for the highest energy bin which was only required to have a significance $\geq 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 (constanttbabs(bbody+powerlawhighecut)) results in $\chi_{\nu}^{2}$=1.05 for 483 dof, while adopting Model 2 (constanttbabspcfabspowerlaw*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 $\chi$ 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$\sigma$ 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 $\chi$ 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 $\chi$ 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$\sigma$ (3.32$\sigma$) above the model prediction while the lowest energy bin of the low state spectrum is 2.72$\sigma$ (3.31$\sigma$) below the model prediction, resulting in a net difference of 6.70$\sigma$ (6.63$\sigma$). If a model is a good fit to the data, we expect the $\chi$ 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 $\chi$ 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 $\chi$ residuals of the two spectra should exceed 4.0$\sigma$ 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 $\chi$ residual differences exceeding 4.0$\sigma$, and the difference of $6.6\sigma-6.7\sigma$ exhibited by the lowest energy bin is much higher than any other bin.

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 $\chi$ 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 $\chi$ residual difference to be $<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^{2}_{\nu}$ 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 $\chi$ 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.

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 $\chi$ 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 §IV.1 are not substantially biased by the fact that the XMM-Newton and NuSTAR observations are not perfectly simultaneous.

V. Discussion

V.1. The physical origin of the soft excess

As discussed in §IV, 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_{\mathrm{H}}\approx 4\times 10^{22}$ cm*-2*, 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_{\mathrm{H}}$ of the partial-covering absorber. Thus, the partial-covering model can naturally explain the observed variations in $N_{\mathrm{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_{\mathrm{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_{\mathrm{H}}$, but found that it can only account for about 25% of the measured variations, which span the range of 3–12$\times 10^{22}$ cm*-2*). 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_{\mathrm{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_{\mathrm{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_{X}<10^{35}$ erg s*-1* (Reig et al. 2009; Bartlett et al. 2013). Furthermore, the flux of the blackbody component of IGR J18214-1318 is 73%$\pm 10$% of the total 0.3–10 keV flux, which is much higher than the $\sim 30$% hot spot blackbody flux fraction typically seen in other HMXBs (e.g. Mukherjee & Paul 2005; La Palombara & Mereghetti 2006; La Palombara & Mereghetti 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$ keV and produces significant emission at $\lesssim 1$ 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.

V.2. The compact object in IGR J18214-1318

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 $<25$ keV argues in favor of an NS since BH HMXBs exhibit power-law spectra out to $\gtrsim 100$ keV (Zdziarski, 2000). Fitting the NuSTAR spectrum above 20 keV with a power-law model, we find that $\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 $\Gamma=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}-10^{13}$ G). The lack of detected pulsations with periods $\lesssim 1$ hrr 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.

V.3. Comparison to other HXMBs

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}}>$1 hr) as a point of comparison.

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$ 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_{\mathrm{H}}\approx 10^{23}$ cm*-2*. 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 $\sim 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_{\mathrm{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_{\mathrm{H}}\sim 10^{24}$ cm*-2*; Walter et al. 2003), which is due to its sgB[e] donor star (Filliatre & Chaty, 2004). The relatively slow 400 km s*-1* 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_{X}\sim 10^{36}$ erg s*-1*) 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$\pm$0.27 $M_{\odot}$ (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 (Reig et al. 2009; Corbet et al. 1999). 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\gtrsim 10^{14}$ G, was slowed down efficiently by the propeller effect before its magnetic field significantly decayed to its current expected value of $\sim 10^{12}$ 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\times 10^{13}$ and $3\times 10^{15}$ 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 (Bonning & Falanga 2005; Torrejón et al. 2004).

Neutron stars 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 $\sim$1.6 hr and $\sim$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 $\sim 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.

VI. Conclusions

Timing analysis of the XMM-Newton and NuSTAR observations of IGR J18214-1318 shows that this HMXB has strong levels of aperiodic variability but no pulsations with periods shorter than an hour. Joint fitting of the XMM-Newton and NuSTAR spectra reveals the presence of an exponential cutoff with e-folding energy $<25$ keV. Thus, although we cannot definitively identify the nature of the compact object in this system, the spectral cutoff energy is a strong indication that it is an NS. A soft excess is also detected in the spectrum of IGR J18214-1318, which may result from the inhomogeneous supergiant wind partially obscuring the X-ray source or blackbody emission with $kT=1.74^{+0.04}_{-0.05}$ keV and $R_{BB}\approx 0.3$ km, which are values typical for emission from hot spots on an NS surface. The partial-covering absorption model is marginally preferred since it can explain the large variations of the total column density measured by different observations; furthermore, the blackbody temperature and the ratio of the blackbody flux to the total X-ray flux measured for IGR J18214-1318 are significantly higher than for other Sg HMXBs, in which blackbody emission seems to arise from the supergiant wind rather than hot spots.

This study demonstrates the usefulness of the combination of XMM-Newton and NuSTAR observations to identify the likely nature of compact objects in HMXBs. Similar observations of other HMXBs will help place better constraints on the BH/NS ratio in HMXB populations. We know of four other HMXBs which, like IGR J18214-1318, have supergiant donor stars, cutoff energies $\lesssim 30$ keV suggesting they harbor NSs, but no detected pulsations despite sensitive timing observations. These sources also resemble 2S 0114+950, an Sg HMXB with 2.7 hr pulsations thought to host a former magnetar. Determining how many of the non-pulsating HMXBs in fact harbor long-period pulsars could shed light on the origins of magnetars and their possible connection to long-period pulsars in HMXBs.

We thank the anonymous referee for suggestions that improved the clarity and thoroughness of this paper. This work made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the NuSTAR Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). In addition, F.M.F. acknowledges support from the National Science Foundation Graduate Research Fellowship. F.M.F. and J.A.T. acknowledge partial support from NASA under XMM Guest Observer grant NNX15AG31G. R.K. acknowledges support from the Russian Science Foundation (grant 14-12-01315). L.N. wishes to acknowledge financial support by ASI/INAF grant I/037/12/0.