NuSTAR DISCOVERY OF A CYCLOTRON LINE IN THE ACCRETING X-RAY PULSAR IGR J16393-4643

During a survey of the Galactic Plane, the ASCA space telescope detected a new X-ray source, AX J163904-4642, in the direction of the Norma spiral arm tangent (Sugizaki et al. 2001). The source was initially classified as a microquasar (Combi et al. 2004) given the absorbed power-law shape of the ASCA-derived spectrum, and given that its X-ray position coincides with counterpart candidates from the radio and infrared bands, as well as with the unidentified gamma-ray source 3EG J1639-4702 (Hartman et al. 1999). Surveying the same region a few years later, the INTEGRAL space telescope detected IGR J16393-4643 which was shown to be the hard X-ray counterpart to the ASCA source (Bird et al. 2004; Malizia et al. 2004).

A follow-up observation of IGR J16393-4643 with XMM-Newton revealed several interesting clues about the source's nature (Bodaghee et al. 2006). First, the refined position from XMM-Newton excluded all previous multi-wavelength counterpart candidates (other than the ASCA source) that were the basis for the microquasar interpretation, and pointed instead to a faint object from the Two Microns All-Sky Survey (Cutri et al. 2003) named 2MASS J16390535-4642137. Thus, the association of IGR J16393-4643 with 3EG J1639-4702 was attributed to a chance alignment of unrelated objects. Second, the spectral properties were typical of a wind-accreting pulsar (e.g., Nagase 1989); i.e., a large absorbing column, a hard X-ray continuum with an exponential cutoff around 20 keV, and iron fluorescence lines. Third, a timing analysis of the XMM-Newton and INTEGRAL data led to the discovery of a coherent pulsation period of 911 s indicative of a slowly rotating, magnetized neutron star. This pulsation period was confirmed in observations taken with RXTE (Thompson et al. 2006), Chandra (Fornasini et al. 2014), and Suzaku (Islam et al. 2015).

These results suggest that IGR J16393-4643 is an obscured high-mass X-ray binary (HMXB) in which a compact object (in this case, a spinning neutron star) accretes the wind shed by a massive donor star (White et al. 1983; Nagase 1989; Bildsten et al. 1997). This view is supported by the subsequent detection of a 4.2 day orbital period in RXTE and Swift data (Corbet et al. 2010; Corbet & Krimm 2013; Coley et al. 2015) which places IGR J16393-4643 within the wind-fed HMXB systems in the pulse-versus-orbital-period diagram of Corbet (1986).

However, the exact spectral class of this donor star is uncertain. The optical/infrared spectrum of 2MASS J16390535-4642137 indicates a spectral class of BIV-V (Chaty et al. 2008), whereas analysis of the K_s-band spectrum of the same object suggests a late-type KM star in a symbiotic binary system (Nespoli et al. 2010). To help clarify this issue, we performed a Chandra observation of IGR J16393-4643 that provided an X-ray position with sub-arcsecond accuracy (Bodaghee et al. 2012). The Chandra position is R.A. $({\rm{J}}2000)={16}^{{\rm{h}}}{39}^{{\rm{m}}}05\overset{{\rm{s}}}{.}47$ and decl. = −46°42'130 with an error radius of 06 (90% confidence). This error circle excludes the 2MASS star and suggests a faint, blended, and likely distant B star that appears in the mid-infrared at wavelengths longer than 5 μm. Indeed, an orbital-period analysis suggests that the donor star is a B giant with a mass greater than 7 M_⊙ (Coley et al. 2015).

In 2014 June, the Nuclear Spectroscope Telescope Array (NuSTAR: Harrison et al. 2013) observed the field of IGR J16393-4643 as part of its survey of the Norma Arm region (Bodaghee et al. 2014; Fornasini et al. 2014). NuSTAR provides exceptional angular (18'' FWHM 58'' half-power diameter) and spectral resolution (400 eV) around 10 keV. In this paper, we present results from these NuSTAR observations of IGR J16393-4643 in addition to a simultaneous snapshot observation taken with Swift . Section 2 describes the analysis of the X-ray data with results from timing and spectral analyses shown in Section 3. Insights into the nature of IGR J16393-4643 are discussed in Section 4.

Table 1 lists details of the observations of IGR J16393-4643 included in this analysis. All data were analyzed using HEASoft 6.16. The NuSTAR data consist of the two focal plane modules A and B (FPMA and FPMB) where each module has a field of view (FOV) of 13' × 13'. Raw event lists from observation ID (ObsID) 30001008002 were reprocessed with nupipeline, which is part of the NuSTAR Data Analysis Software¹³ (NuSTARDAS 1.4.1), while employing the most recent calibration database files available at the time (CALDB: 2014 August 14).

Figure 1 presents the cleaned images for each module in the 3–79 keV energy band. From the cleaned event lists of each module, we extracted source spectra and light curves using a 50''-radius circle centered on the Chandra position for IGR J16393-4643. Around 70% of the source photon energy is enclosed within this radius (Madsen et al. 2015). Given the size of this extraction region, and given that the Chandra position is within 3'' of the brightest pixel in the image from each NuSTAR FPM, we chose not to correct these images for the known systematic offset of coordinates.

Zoom In Zoom Out Reset image size

The bright background feature affecting both modules (Figure 1) is due to unfocused stray-light photons from GX 340+0, an unrelated object situated just outside the FOV. Since the source extraction regions have some fraction of their area contaminated by stray-light photons, we selected background extraction regions (100''-radius circles) that encompassed a similar fraction of stray-light photons. Exposure differences due to vignetting were accounted for in the response matrices and spectra. The effective exposure time at the position of IGR J16393-4643 is 50.579 ks.

During the NuSTAR observation, Swift-XRT (Burrows et al. 2005) also observed the source yielding an effective exposure time of 2.804 ks (ObsID 00080170004). We extracted a spectrum in the 0.5–10 keV energy range that extends the source continuum below the 3 keV limit of NuSTAR thereby enabling the column density to be constrained. The combined Swift-NuSTAR spectra were fit in Xspec 12.8.2 (Arnaud 1996) where we assumed Wilms et al. (2000) abundances and Verner et al. (1996) photo-ionization cross-sections.

3.1. Timing Analysis

Figure 2 presents the source and background light curves binned at 100 s in the 3–79 keV energy range combining counts from NuSTAR FPMA and FPMB. The hardness ratio is defined as $\tfrac{H-S}{H+S}$ where S and H represent count rates in 3–10 keV ("soft") and 10–79 keV ("hard"), respectively. A dividing value of 10 keV allocates a roughly even number of source counts between soft and hard energy bands.

Zoom In Zoom Out Reset image size

The known pulsation is detected at a period of 904.0 ± 0.1 s in the source (+background) light curve with 0.1 s resolution. The best-fitting period is obtained from the fast algorithm for Lomb–Scargle periodograms (Lomb 1976; Scargle 1982) developed by Press & Rybicki (1989), while the error on the pulse period is derived from the analysis methods of Horne & Baliunas (1986) and Leahy (1987) which yield consistent uncertainties. Figure 3 displays the periodogram and Figure 4 shows the phase-folded light curve.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As illustrated in Figure 4, the pulse profile begins with a spike in count rates at phase 0.2–0.3, followed by a dip at phase 0.35, and then a broad secondary peak at phase 0.4–0.85. This bimodal pulse profile is similar to the pulse profiles recorded for this source by INTEGRAL (20–40 keV: Bodaghee et al. 2006) and XMM-Newton (0.3–10 keV: Bodaghee et al. 2006) as well as with RXTE (3–24 keV: Thompson et al. 2006) and Suzaku (0.3–50 keV: Islam et al. 2015).

The pulse fraction $\left(\equiv \tfrac{{I}_{\mathrm{max}}-{I}_{\mathrm{min}}}{{I}_{\mathrm{max}}+{I}_{\mathrm{min}}}\right)$ of 38% ± 1% is consistent with previous measurements by INTEGRAL (57% ± 24%: Bodaghee et al. 2006), XMM-Newton (38% ± 5%: Bodaghee et al. 2006), and Suzaku-XIS (33%: Islam et al. 2015), but it is higher than the fraction measured with RXTE (21% ± 1%: Thompson et al. 2006). The pulse fraction increases with energy (e.g., Lutovinov & Tsygankov 2009) reaching 60% for the NuSTAR 10–79 keV energy band, which is consistent with the pulse fraction from Suzaku-PIN (65% in 12–50 keV: Islam et al. 2015).

In Figure 6, the evolution of the pulsation period is shown as measured by different authors. A least-squares fit to the data results in a slope of $\dot{P}=-2\times {10}^{-8}$ s s⁻¹, which implies that the neutron star is undergoing a long-term spin-up trend at a rate of around −0.6 s per year. This is twice the value of $\dot{P}$ reported by Thompson et al. (2006) using RXTE data. The frequency derivative is $\dot{\nu }=3\times {10}^{-14}$ s⁻². There could be some stochastic variation around the average spin-up value, as seen e.g., in Vela X-1 (Tsunemi 1989; Ikhsanov et al. 2014), but the measurements are too sparse to make a definitive claim.

Zoom In Zoom Out Reset image size

Figure 7 presents the power density spectrum (PDS) of IGR J16393-4643 after subtracting the average pulsed periodic component from the light curve. A broken power-law model applied to the PDS yields a break frequency of 0.00108(21) Hz with Γ₁ = −0.53 ± 0.13 and Γ₂ = −2.08 ± 0.16. The spectral break is necessary given that its inclusion significantly improves the fit quality (χ²/dof = 22.8/23, where dof is the degrees of freedom, compared to χ²/dof = 94.6/25 without the break). In contrast to the recent spin-up detection in 2RXP J130159.6-635806 (Krivonos et al. 2015), where the break in the power spectrum was shifted with respect to the pulsation period, in our case we observe pulsations directly at the break of the power spectrum.

Zoom In Zoom Out Reset image size

A phase lag was suggested by Islam et al. (2015) when they compared the low-energy and high-energy pulse profiles from Suzaku . The minimum of the pulse profile in the low-energy NuSTAR data (3–10 keV) trails the minimum in the high-energy NuSTAR data (10–79 keV) by 0.1 in phase (∼90 s); i.e., consistent with the Suzaku results (Islam et al. 2015). However, unlike in the Suzaku data, the pulse profiles in the NuSTAR data present two maxima and a dip that are synchronized across energy bands (Figure 5).

3.2. Spectral Analysis

The spectral data from the NuSTAR FPMA and FPMB were rebinned so that each energy bin had a significance of at least 4σ. The data were fit simultaneously with an absorbed power law where the instrumental constant was fixed at 1 for FPMA and was allowed to vary for FPMB. In all cases, the instrumental constant that was free to vary remained consistent with 1. The column density was fixed to 2 × 10²³ cm⁻² (Bodaghee et al. 2006). This yielded a poor fit with a reduced ${\chi }_{\nu }^{2}/\mathrm{dof}=7.8/949$ and residuals below 5 keV and above 30 keV. A bremsstrahlung model did not improve the fit by much (${\chi }_{\nu }^{2}/\mathrm{dof}=5.8/949$).

A thermal blackbody (kT = 4.27 ± 0.04 keV) or a power law with an exponential cutoff (Γ = −0.70 ± 0.06 and E_cut = 6.5 ± 0.2 keV) offered better fits to the spectral data (${\chi }_{\nu }^{2}/\mathrm{dof}=1.63/949$ and ${\chi }_{\nu }^{2}/\mathrm{dof}=1.32/948$, respectively). Some residuals remain near 6.4 keV where an iron Kα line is known from XMM-Newton observations (Bodaghee et al. 2006); the addition of a model component for this line leads to a negligible improvement in the fit quality, and so it is not required by the data. The non-detection of this feature is likely due to the lower spectral resolution of NuSTAR at this energy compared with XMM-Newton. Where NuSTAR excels is above 10 keV, and the spectrum shows residuals around 30 keV which suggests a cyclotron resonant scattering feature (CRSF).

To test the significance of this candidate cyclotron line, the NuSTAR spectral data were rebinned so that each energy bin had a significance of least 9σ. We began with the multi-component continuum model from the broadband spectral fits (see below) with the column density fixed to its best-fitting value (N_H = 4 × 10²³ cm⁻²). Adding a cyclotron line to the model reduced the χ² of the fit by 30. In order to estimate the significance of the cyclotron line, we relied on the Bayesian posterior predictive probability value ("ppp-value") as described in Protassov et al. (2002). See Bellm et al. (2014) and Bhalerao et al. (2015) for recent applications of this technique. We determined the reference distribution empirically using a Monte Carlo method where we simulated 1000 trials with the Xspec tool simftest, allowing the centroid energy and width of the cyclotron line to vary within their 90% confidence regions. It is important to note that simftest was used only to simulate the data within a reasonable range of parameter uncertainties and not to perform an F-test. From our simulations, we calculated the change in chi-squared values for the model that includes the cyclotron feature, and for the model without the cyclotron feature (the "null hypothesis"). These simulations returned a maximum change in χ² of 14. The probability of finding the observed change in χ² by chance is 2 × 10⁻⁶ which corresponds to 4.8σ significance (4.3σ after accounting for trials).

Swift spectral data, where each bin contained a minimum of five source counts, were then jointly fit with the NuSTAR data. Using single-component models for the continuum led to fits of insufficient quality and unconstrained spectral parameters. The best fit was obtained with a two-component continuum model incorporating both a soft thermal component (a radial blackbody or bbodyrad) and a hard non-thermal component (a power law with an exponential cutoff or cutoffpl). Photoelectric absorption from molecular hydrogen (tbabs) and a cyclotron absorption line (cyclabs) were included (see also Schönherr et al. 2007, for a more elaborate model).

The fit quality is excellent as attested by the reduced ${\chi }_{\nu }^{2}/\mathrm{dof}$ of 0.99/970 and by the lack of significant residuals as shown in Figure 8. Table 2 lists the spectral parameters of this model. The broadband X-ray flux is about twice as high as was seen previously with Suzaku (Islam et al. 2015) and XMM-Newton (Bodaghee et al. 2006). The effective exposure time is around 51 ks.

Zoom In Zoom Out Reset image size

Table 2.  Spectral Parameters from the Model that Best Fits the Combined Swift-NuSTAR Spectrum of IGR J16393-4643

Model Parameter	Full	Peak	Valley	Unit
const∗tbabs∗cyclabs∗(bbodyrad+cutoffpl)
C	${0.8}_{-0.1}^{+0.2}$	${0.9}_{-0.1}^{+0.2}$	0.5 ± 0.1	⋯
NH	${42}_{-4}^{+3}$	42	42	1022 cm−2
Γ	$-{1.8}_{-0.4}^{+0.6}$	$\leqslant -2.3$	$\sim -2.5$	⋯
${E}_{\mathrm{cut}}$	${5.0}_{-0.4}^{+0.6}$	${4.6}_{-0.2}^{+0.4}$	${4.3}_{-0.3}^{+0.1}$	keV
norm. at 1 keV	3 ± 2	∼1	∼0.4	${10}^{-5}$ ph cm−2 s−1
kT	1.4 ± 0.1	1.6 ± 0.1	1.6 ± 0.1	keV
norm.	0.7 ± 0.2	∼0.6	∼0.6	⋯
${E}_{\mathrm{cyc}}$	${29.3}_{-1.3}^{+1.1}$	29.4 ± 1.2	${30.4}_{-0.9}^{+1.4}$	keV
${\sigma }_{\mathrm{cyc}}$	${4}_{-2}^{+5}$	$\leqslant 5$	${3}_{-1}^{+4}$	keV
depth	0.4 ± 0.1	${0.6}_{-0.3}^{+0.4}$	${1.9}_{-0.6}^{+1.6}$	keV
${F}_{\mathrm{abs}}$	${19.6}_{-1.4}^{+0.3}$	22 ± 11	${12}_{-7}^{+4}$	${10}^{-11}$ erg cm−2 s−1
${F}_{\mathrm{unabs}}$	${28.3}_{-2.0}^{+1.3}$	31 ± 15	18 ± 8	${10}^{-11}$ erg cm−2 s−1
${\chi }_{\nu }^{2}/\mathrm{dof}$	$0.99/970$	$1.06/749$	$1.10/749$	⋯

Note. Spectral parameters are shown for the "full" or time-averaged spectrum, as well as for the 904 s pulse phase-resolved spectra where the "peak" corresponds to counts from phases 0.6–0.85, while the "valley" includes only those counts corresponding to phases in 0.0–1.5 and 0.9–1.0 (see Figure 4). The column density for the phase-resolved spectra is fixed to the value from the "full" spectrum. C is an instrumental cross-calibration coefficient which is fixed at 1 for Swift and variable for NuSTAR . Flux values are given as observed ("abs") and corrected for absorption ("unabs") in the 0.5–80 keV energy range. Uncertainties on the flux offer a more realistic representation of the error range than do the uncertainties on the normalizations which are large and mostly omitted. Errors are quoted at 90% confidence.

Download table as:  ASCIITypeset image

To test for variations in the spectral parameters with respect to the pulse profile, we extracted NuSTAR spectra corresponding to the maxima (or "peak") and minima (or "valley") of the pulse profile. For the peak spectrum, we selected events corresponding to phases 0.6–0.85 (for a count rate of 1.22 ± 0.01 cps in 3–79 keV), whereas events with phases 0–0.15 and 0.9–1 were chosen for the valley spectrum (a count rate of 0.65 ± 0.01 cps). For reference, the "full" or time-averaged source count rate is 1.05 ± 0.01 cps.

Each of these phase-resolved NuSTAR spectra covers one-quarter of the full pulse cycle for a total exposure time of ∼13 ks per spectrum. Each spectrum was rebinned to have a minimum significance of 6σ per bin, and was then jointly fit with the Swift spectrum using the multi-component model described above while holding the column density to its optimal value from the time-averaged spectrum (N_H = 4.2 × 10²³ cm⁻²). A few of the spectral parameters could not be easily constrained and so the following results should be considered with that caveat in mind. The cyclotron line energies from both spectra were statistically consistent with each other: E_cyc = 29.4 ± 1.5 keV for the peak spectrum, and E_cyc = 29.6 ± 1.6 keV for the valley spectrum. The only significant differences were the normalization and flux values which were twice as high in the peak spectrum, i.e., consistent with the difference in count rates. An instrumental cross-calibration constant near 0.5 for the valley spectrum is expected since the NuSTAR count rate is nearly half its average value while the Swift spectrum does not change.

The X-ray band spectrum generated from joint Swift and NuSTAR observations offers the sharpest view yet of the broadband spectral energy distribution (SED) of IGR J16393-4643, especially above ∼10 keV. Thermal blackbody photons with a temperature of 1.41 ± 0.12 keV originating near the surface of the neutron star are re-emitted at higher energies by inverse-Compton scattering off the surrounding electron plasma, and their SED assumes the form a hard power law (${\rm{\Gamma }}=-{1.8}_{-0.4}^{+0.6}$) with an exponential cutoff at ${5.0}_{-0.4}^{+0.6}$ keV. Indeed, a thermal Comptonization model also fits the data providing this physical description of the observed phenomena (e.g., Becker & Wolff 2005).

We discovered a CRSF in the NuSTAR spectrum of IGR J16393-4643 lending further evidence that the compact object in this system is a strongly magnetized neutron star. The centroid energy is ${E}_{\mathrm{cyc}}={29.3}_{-1.3}^{+1.1}$ keV, and it does not change appreciably across the pulse profile (but see, e.g., Coburn et al. 2002). This implies a neutron star magnetic field strength of (2.5 ± 0.1) × 10¹² G, which could be up to 40% larger when gravitational redshift is considered. This E_cyc value is equal to that of X Per and is within the range of observed line energies in accreting X-ray pulsars (e.g., Caballero & Wilms 2012; Walter et al. 2015, and references therein).

The X-ray photons arise at the accretion poles (e.g., Burnard et al. 1991; Becker & Wolff 2005) given that the emission follows a coherent pulsation period seen in the NuSTAR data of IGR J16393-4643, as well as in previous observations. Since the discovery of the pulse period in 2006, each of the five subsequent observations shows IGR J16393-4643 spinning at a higher frequency than during the preceding observation. The change in frequency is $\dot{\nu }=3\times {10}^{-14}$ Hz s⁻¹ which converts to a long-term spin-up trend of −0.6 s per year on a pulsation period of 904.0 ± 0.1 s in the most recent data (this work). Considering that our period derivative is twice the value found by Thompson et al. (2006), and considering that the period measurement of Islam et al. (2015) is not well fit by the least-squares approximation, this indicates that $\dot{\nu }$ is generally increasing but not at a constant rate.

Angular momentum must be imparted to the neutron star to change the spin period in this way (Pringle & Rees 1972; Lamb et al. 1973). The environment of the neutron star in an absorbed HMXB such as IGR J16393-4643 is generally assumed to be a quasi-spherical stellar wind (e.g., Bondi 1952; Walter et al. 2006; Shakura et al. 2012, 2015). Accretion of this wind can, under certain conditions, provide a net change in angular momentum (Shakura et al. 2012). Alternatively, a wind-fed accretion disk can develop around the neutron star in IGR J16393-4643 (Ghosh & Lamb 1979; Boerner et al. 1987), and since a BIII-type donor star is not expected to fill its Roche lobe, this favors a transient scenario for the disk. Transient accretion disks have been suggested in other wind-fed pulsars that exhibit long-term spin-up trends (e.g., GX 301-2: Koh et al. 1997).

Naturally, the strong magnetic field of the pulsar influences the dynamics of the accretion disk out to a certain distance (Ghosh & Lamb 1979). Material within the corotation radius is forced to corotate with the neutron star and if its angular velocity is less than the local Keplerian velocity, then the material can be accreted along the field lines (this implies a magnetospheric radius smaller than the corotation radius). Depending on whether the accreted material's angular momentum has the same or an opposite direction as the neutron star's spin, this will cause the neutron star to spin faster or slower, respectively (Waters & van Kerkwijk 1989). The consistent downward trend of the pulsation period in IGR J16393-4643 suggests that material in the accretion disk orbits in the same direction as the neutron star's spin.

Additional accretion dynamics can be elucidated from the PDS. Going from low frequency to high frequency, the PDS changes from a flat power law (slope ∼ −0.5) to a steeper power law (slope ∼ −2). The PDS breaks at a frequency of 0.00108(21) Hz which is very close to the spin frequency of the neutron star (0.0011062(1) Hz). These characteristics of the PDS in persistent, accreting X-ray pulsars are described by Revnivtsev et al. (2009), and they suggest that the neutron star in IGR J16393-4643 is spinning close to corotation with the inner edge of an accretion disk that has been truncated by the magnetic field. Using canonical values for the neutron star mass ($1.4\;{M}_{\odot }$) and magnetic dipole moment (μ = 10³⁰ G cm⁻³ for B ∼ 10¹² G), this gives estimates of 1.1 × 10⁸ cm and 2.0 × 10⁸ cm for the magnetospheric radius assuming mass accretion rates of 10⁻⁶ M_⊙ yr⁻¹ and 10⁻⁷ M_⊙ yr⁻¹, respectively. However, an alternative model for quasi-spherical wind accretion developed by Shakura et al. (2012) could also produce the observed change in spin frequency, and in this case, an accretion disk is not required.

A NuSTAR observation of IGR J16393-4643 has revealed valuable insights into the nature of this source. The detection of a CRSF at ${29.3}_{-1.3}^{+1.1}$ keV allowed us to constrain the magnetic field strength to B = (2.5 ± 0.1) × 10¹² G. This is the first time that the magnetic field has been measured in this object. The cyclotron line was not detected in previous excursions into this energy range by RXTE-ASM, INTEGRAL-ISGRI, and Suzaku-HXD. This result is a testament to the spectral sensitivity and resolving power of NuSTAR.

The pulsation period of the neutron star in IGR J16393-4643 is now at 904.0 ± 0.1 s. Looking at the five measurements made between 2006 and 2014, we find that this period has consistently gotten shorter with time at an average long-term spin-up rate of −0.6 s per year ($\dot{\nu }=3\times {10}^{-14}$ Hz s⁻¹). The slope of the PDS breaks near the pulsation frequency as expected for persistently emitting accretion-powered pulsars. This could indicate that a transient and magnetically truncated accretion disk is almost in corotation with the neutron star whose magnetospheric radius is around 2 × 10⁸ cm, although accretion from a quasi-spherical wind could also lead to the observed change in pulsation frequency.

IGR J16393-4643 is an archetype of the class of heavily obscured wind-accreting pulsars, with its ks-long pulsation period and its column density that is at least an order of magnitude greater than the expected line of sight value. Yet its source classification has been subject to multiple changes over the years and the identity of the donor star remains elusive. The results we obtained thanks to NuSTAR represent new puzzle pieces to add to the still-developing picture we have of this intriguing source.

The authors are grateful to the anonymous referee for their constructive criticism that helped improve the quality of the manuscript. A.B. thanks Dr. Konstantin Postnov. F.M.F. acknowledges support from the National Science Foundation Graduate Research Fellowship. R.K. acknowledges support from from Russian Science Foundation (grant 14-22-00271). The scientific results reported in this article are based on 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; data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC) provided by NASA's Goddard Space Flight Center; NASA's Astrophysics Data System Bibliographic Services; and the SIMBAD database operated at CDS, Strasbourg, France.

Facilities: NuSTAR - The NuSTAR (Nuclear Spectroscopic Telescope Array) mission, Swift - Swift Gamma-Ray Burst Mission.