The X-Ray Pulsar XTE J1858+034 Observed with NuSTAR and Fermi/GBM: Spectral and Timing Characterization plus a Cyclotron Line

As observed by NuSTAR in 2019 November, XTE J1858+034 clearly shows a hard spectrum. The FPMA and FPMB spectra have been fitted simultaneously, allowing for a cross-normalization factor. Although the cross-normalization factor between FPMA and FPMB is usually of the order of a few percent (Madsen et al. 2015), the limited ARF extraction region adopted in our analysis for FPMB (see Section 2) is expected to reduce the cross-normalization value significantly. For the spectral fit, standard phenomenological and semiempirical continuum models have been employed, namely, two variants of the cutoff power-law model (cutoffpl and highecut ∗ pow in XSPEC) and a Comptonization model of soft photons in a hot plasma (compTT in XSPEC; Titarchuk 1994), respectively. To obtain an acceptable fit, the cutoffpl and highecut ∗ pow models need an additional component in the lower energy band, which has been modeled as a blackbody emission as found in other accreting XRPs (see, e.g., La Palombara & Mereghetti 2006). However, the blackbody temperature is high with respect to other XRPs, indicating that the phenomenological model is likely inadequate. Moreover, we also tested a purely physical model of thermal and bulk Comptonization of the seed photons produced by cyclotron cooling (Ferrigno et al. 2009; bwcycl in XSPEC). For a fixed value of mass and radius of the accreting NS, the bwcycl model has six free parameters, namely, the accretion rate $\dot{M}$, magnetic field strength B, accretion column radius r₀, electron temperature T_e , photon diffusion parameter ξ, and Comptonization parameter δ. This model was successfully used to fit the broadband energy spectrum of a number of bright (≳10³⁷ erg s⁻¹) accreting XRPs (see, e.g., Wolff et al. 2016; D'Aì et al. 2017; Epili et al. 2017).

For all tested models, the photoelectric absorption component and elemental abundances were set according to Wilms et al. (2000; tbabs in XSPEC) to account for photoelectric absorption by neutral interstellar matter (or column density N_H) and assuming model-relative (wilm) solar abundances. Given that the Galactic N_H in the direction of the source is about 1.7 × 10²² cm⁻² (HI4PI Collaboration et al. 2016), all models show important local absorption values. All tested models also were equipped with a Gaussian emission line at 6.4 keV to account for the Fe Kα fluorescence emission.

All fit continuum models show absorption-like residuals in the range 40–50 keV. These residuals can be modeled with a Gaussian absorption line (see Figure 3). The improvement in the best-fit statistics is maximum in the compTT model, i.e., Δχ² = 878. Other models show an improvement of Δχ² ⪆ 400, with the lowest Δχ² derived from the highecut∗pow model. The significance of the line in the compTT model has been assessed through Monte Carlo simulations. For this task, the XSPEC simftest routine was adopted, which allows one to simulate a chosen number of spectra based on the actual data and test the resulting Δχ² between each fit when the additional model component (the Gaussian absorption line, in our case) is included. Following Bhalerao et al. (2015) and Bodaghee et al. (2016), the column density parameter was fixed to its best-fit value, and the energy and width of the Gaussian absorption line were left free to vary within their 90% confidence region in order to improve the speed and convergence of the fits. Simulation results are reported in Figure 4 for a 10⁴ iteration process and confirm the significance of the absorption feature at >3σ c.l. Following MarcuCheatham et al. (2015), we also investigated the impact of a variable background normalization on the absorption feature parameters. Using the XSPEC tool recorn, it was found that the absorption line parameters do not change significantly if the normalization of the background spectrum is increased by up to 50%, thus strengthening the interpretation of the absorption feature as real and not due to artifacts. The feature was also observed in phase-resolved spectra presented in the accompanying paper by Tsygankov et al. (2021). Interpreting the feature at 48 keV as a CRSF and assuming a gravitational redshift of z_g = 0.3 (for an NS mass and radius of 1.4 M_⊙ and 10 km, respectively), a magnetic field strength of B = (5.4 ± 0.1) × 10¹² G is obtained.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Following the bwcycl model instructions, ¹⁴ it is convenient to freeze some of the model parameters in order to improve the computational speed and help the fit converge to the best-fit parameters. Once the best fit was found, the column density N_H was also fixed to its best-fit value to help the fit converge and obtain parameter errors. The mass and radius of the NS were fixed to their canonical values of 1.4 M_⊙ and 10 km, respectively. However, it is preferable to also fix the values of the NS magnetic field, its distance, and its mass accretion rate (as derived by the observed luminosity). For XTE J1858+034, there are no previous conclusive estimations of the magnetic field, while a measurement of the distance is necessary for the latter two parameters. As mentioned in Section 1, the closest Gaia counterpart to the nominal X-ray position found by Molkov et al. (2004) was unlikely to be associated with the X-ray source or the optical counterpart proposed by Reig et al. (2005). This was ascertained by Tsygankov et al. (2021), who showed that either of those possible counterparts is consistent with the much better constrained X-ray source location available through new Chandra observations. However, the distance to the X-ray system had not been estimated before their work, which, in any case, did not constrain it very much. Therefore, we opted for a different approach to obtain a more stringent value of the distance, based on the spin-up ($\dot{P}$) measured by the Fermi Gamma-ray Burst Monitor (GBM).

To this aim, the publicly available spin-frequency values from the GBM were used. ¹⁵ The spin-up was measured during an interval of about 6 days around MJD 58,786. The resulting spin-up value is $| \dot{P}| =10.5421(6)$ s yr⁻¹ (see also Malacaria et al. 2020 and references therein). Since the orbital parameters of this system are unknown, we tested the contribution of orbital modulation to the observed spin-up. First, although Doroshenko et al. (2008) reported a possible orbital period for this source of about 380 days, we notice that its significance is low, while a visual inspection of the Swift/BAT (Krimm et al. 2013) data ¹⁶ for this source revealed an outburst recurrence of ∼81 days (see Figure 5), here assumed as the orbital period. Moreover, for a K- or M-type optical companion star (Tsygankov et al. 2021), we adopted a value of the mass function $f{(M)=M{{}_{* }}^{3}{\sin }^{3}i/({M}_{* }+{M}_{\mathrm{NS}})}^{2}\,=\,1$, where M_* and M_NS are the mass of the companion and NS, respectively, and i is the binary system inclination. This corresponds to a value of the semimajor projected axis, ${a}_{x}\,\sin \,i\simeq 190\,$ lt-s. Also, the data required only a small value of the eccentricity, assumed here as e ≃ 0.1. Finally, an epoch of T₀ = 53,436 MJD was chosen at the beginning of the first of the recurring outbursts. An argument of periapse of ω = 107° was found to best fit the GBM frequency values assuming no accretion torque. This orbit assumed, the maximum orbital contribution to the spin-up was found to be only about 10%.

Zoom In Zoom Out Reset image size

Assuming a magnetic field strength of 5.4 × 10¹² G and adopting the NuSTAR measured flux of 1.5 × 10⁻⁹ erg cm⁻² s⁻¹ (see Tables 1 and 2), the standard accretion-disk torque theory (Ghosh & Lamb 1979) can be used to infer the distance of the source according to the equation

$$\begin{eqnarray}\begin{array}{rcl}-\dot{P} & = & 5\times {10}^{-5}\,{\rm{s}}\,{\mathrm{yr}}^{-1}\,{M}_{1.4}^{-3/7}\,{R}_{6}^{6/7}\,{I}_{45}^{-1}\\ & & \times \,n({\omega }_{s}){\mu }_{30}^{2/7}{\left(P{L}_{37}^{3/7}\right)}^{2},\end{array}\end{eqnarray} \tag{ 1 }$$

where M_1.4 is the NS mass in units of 1.4 M_⊙, R₆ is the NS radius in units of 10⁶ cm, I₄₅ is the moment of inertia in units of 10⁴⁵ g cm², μ₃₀ is the magnetic moment in units of 10³⁰ G cm³, n(ω_s ) ≈ 1.4 is the dimensionless torque, P is the spin period in seconds, and L₃₇ is the bolometric luminosity in units of 10³⁷ erg s⁻¹. For a measured $| \dot{P}| =10.5\,$ s yr⁻¹, Equation (1) allows one to infer a distance of d = 10.9 ± 1.0 kpc (estimated uncertainty at 1σ c.l.). This distance value is also independently confirmed by the analysis of the optical companion star, as reported in the accompanying paper by Tsygankov et al. (2021), and it was used to characterize different configurations of the bwcyc model. The corresponding mass accretion rate $\dot{M}=1.2\times {10}^{17}\,$ g s⁻¹ was adopted altogether, derived assuming a luminosity $L=\eta \dot{M}{c}^{2}$ with efficiency η = 0.2 (Sibgatullin & Sunyaev 2000). The different configurations of the tested model are reported in Table 2, including a set with the magnetic field strength as a free parameter (BWCYCa), one with the magnetic field strength fixed to 5.4 × 10¹² G (BWCYCb), and one with a free $\dot{M}$ and fixed magnetic field strength (BWCYCc).

Table 1. Best-fit Results of XTE J1858+034 Spectral Analysis with a Cutoff Power-law Model cutoffpl and a Comptonization Model compTT

	cutoffpl	compTT
CFPMB	${0.747}_{-0.001}^{+0.001}$	${0.747}_{-0.001}^{+0.001}$
NH (1022 cm−2)	${7.6}_{-0.4}^{+0.4}$	${5.8}_{-0.3}^{+0.3}$
kTbb (keV)	${5.2}_{-0.1}^{+0.2}$	⋯
Normbb	${0.0152}_{-0.0007}^{+0.0002}$	⋯
EKα (keV)	${6.47}_{-0.02}^{+0.02}$	${6.48}_{-0.02}^{+0.02}$
σKα (keV)	${0.26}_{-0.02}^{+0.02}$	${0.28}_{-0.02}^{+0.02}$
NormKα (10−4)	${5.5}_{-0.4}^{+0.4}$	${5.9}_{-0.4}^{+0.4}$
Γ	${0.03}_{-0.22}^{+0.29}$	⋯
highecut (keV)	${3.5}_{-0.5}^{+1.5}$	⋯
NormΓ a	${0.025}_{-0.004}^{+0.006}$	⋯
T0 (keV)	⋯	${1.02}_{-0.02}^{+0.02}$
kTcompTT (keV)	⋯	${5.61}_{-0.04}^{+0.05}$
τp (keV)	⋯	${7.07}_{-0.06}^{+0.06}$
NormcompTT	⋯	${0.0228}_{-0.0003}^{+0.0003}$
Egabs (keV)	${46.8}_{-0.8}^{+1.0}$	${48.0}_{-0.7}^{+0.8}$
σgabs (keV)	${7.7}_{-0.7}^{+1.0}$	${8.6}_{-0.5}^{+0.6}$
Strength gabs	${14.9}_{-2.6}^{+5.1}$	${21.3}_{-2.5}^{+2.0}$
Flux b	${1.499}_{-0.003}^{+0.003}$	${1.499}_{-0.003}^{+0.003}$
χ2/dof	1645/1574	1643/1573

Notes. All reported errors are at 1σ c.l.

^a In units of photons⁻¹ keV⁻¹ cm⁻² s⁻¹ at 1 keV. ^b Flux calculated for the entire model in the 3–60 keV band and reported in units of 10⁻⁹ erg cm⁻² s⁻¹. Flux values with estimated errors were derived using the cflux model from XSPEC as resulting from FPMA.

Download table as:  ASCIITypeset image

Table 2. Best-fit Results of XTE J1858+034 Spectral Analysis with Different Configurations of the Physical Bulk+Thermal Comptonization Model bwcyc

	BWCYCa	BWCYCb	BWCYCc
CFPMB	${0.748}_{-0.001}^{+0.001}$	${0.748}_{-0.001}^{+0.001}$	${0.748}_{-0.001}^{+0.001}$
NH (1022 cm−2)	${8.6}_{-0.2}^{+0.2}$	${8.4}_{-0.2}^{+0.2}$	${8.4}_{-0.8}^{+0.3}$
EKα (keV)	${6.47}_{-0.02}^{+0.02}$	${6.48}_{-0.02}^{+0.02}$	${6.47}_{-0.02}^{+0.02}$
σKα (keV)	${0.27}_{-0.02}^{+0.02}$	${0.28}_{-0.02}^{+0.02}$	${0.27}_{-0.02}^{+0.02}$
NormKα (10−4)	${5.5}_{-0.4}^{+0.4}$	${5.6}_{-0.3}^{+0.3}$	${5.6}_{-0.4}^{+0.4}$
ξ	${4.2}_{-1.7}^{+0.7}$	${3.2}_{-0.8}^{+0.8}$	${2.9}_{-0.6}^{+1.0}$
δ	${0.4}_{-0.1}^{+0.3}$	${0.6}_{-0.1}^{+0.3}$	${0.7}_{-0.2}^{+0.3}$
B (1012 G)	${4.4}_{-0.3}^{+0.5}$	5.4 (fixed)	5.4 (fixed)
$\dot{M}$ (1017 g s−1)	1.2 (fixed)	1.2 (fixed)	${1.1}_{-0.1}^{+0.1}$
Te (keV)	${5.6}_{-0.6}^{+0.2}$	${5.5}_{-0.5}^{+0.1}$	${5.2}_{-0.4}^{+0.2}$
r0 (m)	${73}_{-22}^{+9}$	${66}_{-12}^{+22}$	${49}_{-10}^{+29}$
d (kpc)	10.9 (fixed)	10.9 (fixed)	10.9 (fixed)
Egabs (keV)	${48.3}_{-0.7}^{+0.7}$	${48.6}_{-1.3}^{+0.5}$	${48.3}_{-1.1}^{+1.0}$
σgabs (keV)	${10.3}_{-2.1}^{+0.7}$	${9.6}_{-1.2}^{+1.2}$	${9.3}_{-0.6}^{+1.3}$
Strengthgabs	${25.4}_{-6.4}^{+4.2}$	${27.6}_{-4.2}^{+3.1}$	${25.3}_{-3.9}^{+4.0}$
Flux b	${1.499}_{-0.003}^{+0.003}$	${1.499}_{-0.002}^{+0.003}$	${1.499}_{-0.003}^{+0.003}$
χ2/dof	1648/1574	1648/1575	1648/1574

Notes. All reported errors are at 1σ c.l.

^a Unconstrained. ^b Flux calculated in the 3–60 keV band and reported in units of 10⁻⁹ erg cm⁻² s⁻¹. Flux values with estimated errors were derived using the cflux model from XSPEC as calculated for FPMA.

Download table as:  ASCIITypeset image

For the timing analysis, the nuproducts task was used to obtain light curves out of calibrated and cleaned events. These light curves were corrected for live-time, exposure, and vignetting effects and extracted in the following energy bands: 3–10, 10–20, 20–30, 30–40, 40–60, and 3–60 keV.

All light curves were barycentered using the barycorr tool and the NuSTAR clock correction file nuCclock20100101v103. The light curve in the 3–60 keV energy band was binned to 5 s and used to search for pulsations around the known 221 s periodicity with the epoch folding method (Leahy 1987). The procedure results in a measured period of P = 218.393(2) s. Pulsations are significant at >99%. The uncertainty was estimated by simulating 500 light curves based on real data and altered with Poisson noise.

Light curves in different energy bands were folded to the best-fit spin period to obtain pulse profiles with a resolution of 20 phase bins (see Figure 6). In turn, these were used to explore the pulsed fraction variation as a function of the energy (see Figure 7). The pulsed fraction here is defined as $({I}_{\max }-{I}_{\min })/({I}_{\max }+{I}_{\min })$, where ${I}_{\max },{I}_{\min }$ are the maximum and minimum pulse profile count rate, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The hard spectrum of XTE J1858+034 resembles that of other accreting XRPs observed at both low and high luminosity and well fit by a compTT model (Mukerjee et al. 2020 and references therein). The observation of compTT spectra in accreting XRPs is usually interpreted as the result of thermal Comptonization processes in which the thermal energy of the accreting gas is transferred to the seed photons originating from the NS hot spots (Becker & Wolff 2007). An increasing number of these sources also show that an additional compTT component emerges in the high-energy range of the spectrum at low-luminosity stages (Tsygankov et al. 2019a, 2019b). Although the formation of such a component is not clear yet, it is likely due to a combination of cyclotron emission and following thermal Comptonized emission from a thin overheated layer of the NS atmosphere (see Tsygankov et al. 2019b, and references therein). In this context, X Persei is a remarkable case, since it has been shown that the cyclotron line in its spectrum can be mimicked by the convolution of the two compTT spectral components around the energy where the flux from the low- and high-energy components is comparable (Doroshenko et al. 2012). However, among the sources whose spectrum is formed by two compTT components, X Persei is the one with the highest electron temperature of the hard-energy compTT component, kT ∼ 15 keV. If the absorption feature at ∼48 keV in XTE J1858+034 does in fact result from the blend of two compTT components, the high-energy compTT would peak around 22 keV. This would make XTE J1858+034 the most extreme among the XRPs that show such a spectral shape. However, such a spectral shape has so far only been observed in low-luminosity XRPs, while our analysis (see Section 3.1) shows that the source is located at a relatively large distance of 10.9 kpc, thus implying a high-luminosity source (also supported by the analysis of Tsygankov et al. 2021). A second Gaussian or compTT component that peaked above 45 keV was also tested in place of the absorption feature but could not be successfully fit (${\chi }_{\mathrm{red}}^{2}\gt 1.4$), possibly due to the lack of statistics above 60 keV.

With the newly inferred source distance value of d = 10.9 kpc, the derived flux of 1.5 × 10⁻⁹ erg cm⁻² s⁻¹ implies a luminosity of 2.1 × 10³⁷ erg s⁻¹. Adopting this distance value, all different sets of the bwcyc model are statistically equivalent, and all parameters show acceptable values. For the BWCYCa configuration, the returned magnetic field strength of 4.2 × 10¹² G is consistent within 2σ with that inferred from the candidate CRSF. Notably, the BWCYCb model with fixed values of the distance, magnetic field strength, and accretion rate also fits the data and returns acceptable values of the best-fit parameters. The BWCYCc fit returns the smallest (but still acceptable) r₀ value and a mass accretion rate $\dot{M}$ that is almost coincident with that inferred from the X-ray (isotropic) luminosity.

In any case, when interpreting the results from the bwcyc model, it is important to keep in mind that, as reported in Ferrigno et al. (2009), the BWcyc model may need adjustments in the spectral parameters with respect to the original prescriptions. For example, the best-fit magnetic field value may differ from that inferred by the CRSF if the spectrum is formed at an NS site that is spatially different from the CRSF-forming region. Likewise, the best-fit mass accretion rate $\dot{M}$ can be different from that inferred by the X-ray luminosity due to an uncertain efficiency conversion factor (η) and anisotropic emission.

The pulse profiles of XTE J1858+034 as observed by NuSTAR show a single-peak structure and a shape that is only weakly energy-dependent (see Figure 6). This is typically observed at low mass accretion rates (see, e.g., Malacaria et al. 2015) and qualitatively interpreted as the beaming pattern resulting from a pencil-beam emission. However, single-peaked pulse profiles are also observed at high accretion rates, like in the case of pulsating ultraluminous X-ray sources, e.g., Swift J0243.6+6124, where single-peak pulse profiles persist at high luminosity and only switch to more complex profiles at super-Eddington luminosity (Wilson-Hodge et al. 2018).

The pulsed fraction shows a considerable energy dependence and almost doubles from 20% in the 3–10 keV energy band to about 40% in the 30–40 keV energy band, above which the lack of statistics prevents us from drawing firm conclusions.