SPECTRAL CHANGES IN THE HYPERLUMINOUS PULSAR IN NGC 5907 AS A FUNCTION OF SUPER-ORBITAL PHASE

NuSTAR (Harrison et al. 2013) data were extracted using the standard nustardas pipeline v1.6.0 as distributed with HEASOFT V6.19. We applied the CALDB 20160706. Source data were extracted from 50'' regions centered on the source and within the pointing uncertainty for the J2000 coordinates of NGC 5907 ULX1. Background spectra were extracted from annuli centered on the same coordinates with inner radii of 90'' and outer radii of 200''. The source becomes background dominated around 10 keV and 20 keV for the 2013 and 2014 data, respectively. We therefore carefully checked that our choice of background region does not influence the results, finding consistent results when using a circular background region with a radius of 100'' located elsewhere on the same chip. Note that we only use the second observation in 2013 because the source was not detected in the first one (W15).

We additionally extracted data taken in SCIENCE_SC (mode 06), during which the optical bench star tracker does not provide a solution for aspect reconstruction. While the source is therefore not reconstructed accurately on the sky, the responses can still be calculated correctly. We used regions of the same size centered on the centroid coordinates of the visible point source in these data sets. See Fürst et al. (2016a) and Walton et al. (2016c) for details on mode 06 extraction. Using these data, we increased the effective exposure time by 10% and 15% for the 2013 and 2014 epoch, respectively. This extra exposure time is included in Table 1.

The 2014 data were obtained during two separate observations, about three days apart (Table 1). We initially analyzed each observation separately but did not find any significant differences between them. We therefore added both observations and treat them as one epoch for the remainder of this paper.

The XMM-Newton (Jansen et al. 2001) observations were reduced with the XMM-Newton Science Analysis System (SAS) v14.0.0 following standard procedures. The raw data files were filtered using epchain and emchain to produce cleaned event lists for each of the EPIC-pn and EPIC-MOS detectors, respectively (Strüder et al. 2001; Turner et al. 2001). As recommended, we use only single and double events for EPIC-pn, and single to quadruple events for MOS. We exclude periods of high background flaring. Science products were then produced using xmmselect, with the source emission extracted from circular regions of radius ∼20''–30'' (depending on the source brightness and its proximity to bad detector columns) and the background estimated from larger areas on the same CCD free of other contaminating point sources. Redistribution matrices and auxiliary response files were generated with rmfgen and arfgen, respectively. After performing the data reduction separately for each of the two MOS detectors, we combined the data from these detectors into a single spectrum for each epoch using the FTOOL addascaspec.

Chandra (Weisskopf et al. 2000) observed NGC 5907 ULX1 for two back-to-back observations in 2012 (Table 1). We extracted the ACIS-S with the standard CIAO v4.8 pipeline. The source spectra were extracted from circular regions with 3'' radius centered on the J2000 coordinates, the background from 25'' radius regions to the northeast from a source-free area. The spectra of the two observations were added with the CIAO tool combine_spectra. In addition, the Chandra observations were performed within a few days of the 2012 XMM-Newton observations and did not show significant changes in spectral shape or flux. We therefore treat all of the 2012 data as a single epoch in the following analysis.

The 2014 data clearly show a hard spectrum with a visible turnover at high energies (Figure 2). Neither a pure power-law (${\chi }^{2}=1367/667$ dof) nor a simple multicolor blackbody of a geometrically thin, optically thick disk accretion disk (diskbb, ${\chi }^{2}=997/667$ dof) adequately describe this shape. This is very similar to other ULXs studied by NuSTAR and to previous studies of NGC 5907 ULX1 (Sutton et al. 2013; Israel et al. 2016a).

Zoom In Zoom Out Reset image size

A phenomenological cutoff power-law model (cutoffpl), however, provides an acceptable fit with ${\chi }^{2}=739$ for 666 dof. This model can be improved by adding a diskbb model with a temperature of ∼0.3 keV (${\chi }^{2}=723$ for 664 dof, ${\rm{\Delta }}{\chi }^{2}=16$ for two additional parameters). Such a model is often used to describe Galactic binaries (e.g., McClintock & Remillard 2006). However, while the NGC 5907 ULX1 spectrum is dominated by the power-law component, similar to the low/hard state of black hole binaries, the values we obtain are very different, e.g., the photon index is much harder (${\rm{\Gamma }}={0.83}_{-0.15}^{+0.13}$ instead of 1.4–1.8) and the folding energy is much lower (${E}_{\mathrm{fold}}={5.3}_{-0.6}^{+0.7}$ keV instead of $\gg 20\,\mathrm{keV};$ see, e.g., Fürst et al. 2015).

A more physically motivated model is the diskpbb model (Mineshige et al. 1994), which allows for a variation of the temperature gradient p of the multicolor blackbody, with $T(r)\propto {r}^{-p}$. Recent NuSTAR results have shown that ULX spectra are often well described with a temperature gradient somewhat shallower than the canonical p = 0.75 expected for a thin disk (e.g., Bachetti et al. 2013; Brightman et al. 2016). Such shallower gradients are expected in sources with very high to super-Eddington luminosities, in which the accretion disk can increase its geometrical thickness due to radiation pressure and advection (Abramowicz et al. 1988). We obtain a comparable fit to the cutoff power-law model with ${\chi }^{2}=740$ for 666 dof. We find $p={0.598}_{-0.010}^{+0.011}$ and ${T}_{\mathrm{in}}={3.46}_{-0.15}^{+0.17}$ keV.

However, this model leaves significant residuals at the highest energies. Such a hard energy excess is expected if a significant fraction of the thermal photons are Compton scattered to higher energies, resulting in an additional high-energy power-law continuum. Evidence for an additional high-energy power-law continuum has now been observed in several ULXs (e.g., Walton et al. 2013, 2014, 2015b; Mukherjee et al. 2015). Here we model the hard excess with the simpl model (Steiner et al. 2009), which emulates up-scattering of the thermal seed photons into a power-law tail. We find an excellent fit with ${\chi }^{2}=725$ for 664 dof (${\rm{\Delta }}{\chi }^{2}=15$ for two additional parameters, Table 2). We find a scattering fraction (i.e., flux in the power-law tail) of ${F}_{\mathrm{sctr}}\approx 0.09;$ however, this value is highly degenerate with the photon index. The corresponding ${\chi }^{2}$-landscape is complex so that a simple uncertainty estimation is not possible. By using the XSPEC steppar command in the 2D-space between ${F}_{\mathrm{sctr}}$ and Γ, we find that ${F}_{\mathrm{sctr}}$ has a lower limit of 0.036 at 90% confidence. We do not find an upper limit due to the fact that Γ can become very high, i.e., steep.

Table 2.  Best-fit Model Parameters for All Epochs

	2014		2013	2012	2003
Parameter	Cpl+Diskbb	Simpl(Diskpbb)	Simpl(Diskpbb)	Diskpbb	Diskpbb
${N}_{{\rm{H}}}\,({10}^{22}\,{\mathrm{cm}}^{-2})$	${0.85}_{-0.12}^{+0.14}$	0.70 ± 0.04	0.65 ± 0.08	0.57 ± 0.10	0.76 ± 0.06
${ \mathcal F }\,({10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1})$ a	${2.54}_{-0.17}^{+0.26}$	${2.44}_{-0.32}^{+0.08}$	${0.86}_{-0.24}^{+0.09}$	0.83 ± 0.05	1.93 ± 0.08
${A}_{\mathrm{disk}}$ b	${3.0}_{-2.3}^{+7.1}$	$({8.9}_{-2.8}^{+3.8})\times {10}^{-4}$	$({2.8}_{-1.6}^{+6.5})\times {10}^{-4}$	$({2.6}_{-2.0}^{+4.5})\times {10}^{-4}$	$({4.9}_{-2.9}^{+4.6})\times {10}^{-4}$
Γ	${0.83}_{-0.15}^{+0.13}$	${1.1}_{-0.0}^{+1.7}$	${1.1}_{-0.0}^{+1.9}$	⋯	⋯
${E}_{\mathrm{fold}}\,(\mathrm{keV})$	${5.3}_{-0.6}^{+0.7}$	⋯	⋯	⋯	⋯
${T}_{\mathrm{in}}\,(\mathrm{keV})$	${0.30}_{-0.06}^{+0.09}$	${2.94}_{-0.30}^{+0.23}$	${2.7}_{-0.9}^{+0.5}$	${3.5}_{-0.7}^{+1.3}$	${3.1}_{-0.4}^{+0.7}$
p	⋯	${0.610}_{-0.013}^{+0.015}$	${0.552}_{-0.020}^{+0.028}$	${0.68}_{-0.05}^{+0.07}$	${0.585}_{-0.020}^{+0.022}$
${F}_{\mathrm{sctr}}$	⋯	$\gt 0.036$	$\lt 0.11$ c	⋯	⋯
${{ \mathcal L }}_{17}\,({10}^{40}$ erg s−1)d	${8.8}_{-0.6}^{+0.9}$	${8.51}_{-1.12}^{+0.28}$	${3.00}_{-0.82}^{+0.30}$	${2.90}_{-0.17}^{+0.18}$	${6.72}_{-0.26}^{+0.28}$
${\chi }^{2}/\mathrm{dof}$	723.52/664	725.49/664	268.09/231	197.78/214	382.06/411
${\chi }_{\mathrm{red}}^{2}$	1.090	1.093	1.161	0.924	0.930

Notes.

^aFlux between 0.5 and 10 keV. ^bNormalization of disk model in units of ${({R}_{\mathrm{in},\mathrm{km}}/{d}_{10})}^{2}\cos (\theta )$. ^cFor a fixed value of ${\rm{\Gamma }}=1.1$. ^dLuminosity between 0.5 and 10 keV for a distance of 17.1 Mpc.

Download table as:  ASCIITypeset image

In all of these models, a small excess around 1.5 keV is visible, which can be linked to larger calibration uncertainties around the known "silicon bump" (Read et al. 2014).

To put the 2014 data into context with previous observations of NGC 5907 ULX1, we perform our own analysis of the 2003, 2012, and 2013 epochs (Figure 1, Table 1). We fit the 2013 broadband data with the same simpl × diskpbb model used for the 2014 data. For the XMM-Newton-only observations in 2003 and the XMM-Newton plus Chandra observation in 2012, we did not include the simpl model because the lack of coverage at high energies does not allow us to constrain its parameters. As can be seen in Figure 2, the hard tail modeled by the simpl model only becomes relevant above ∼15 keV. We also did not use the diskbb + cutoffpl model, as the power-law parameters could only be very weakly constrained with the soft data alone.

The best-fit values for all epochs are given in Table 2. The 2012 and 2013 observations have a very similar 0.5–10 keV flux, around $8.4\times {10}^{-13}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, while both the 2003 and 2014 observations have a higher flux ($\sim 2.5\times {10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$). However, neither the disk temperature nor the temperature gradient p show a clear correlation with flux. Instead, the 2012 data show a significantly higher value of p than the other observations, which are all consistent with each other. With respect to the temperature, all observations are consistent within their uncertainties.

The difference can also be seen in Figure 3, where we show the unfolded spectra of all four epochs. While the observations in 2012 and 2013 show very similar fluxes, their spectral shapes are distinctly different, with the 2012 flux rising much more steeply with energy. On the other hand, the 2003 and 2014 data agree very well with each other.

Zoom In Zoom Out Reset image size

We find a tentative trend of the spectral parameters as a function of the 78-day phase. We list the phases of each epoch in Table 1 and plot the values of ${T}_{\mathrm{in}}$ and p versus flux and phase, respectively, in Figure 4. The 2003, 2013, and 2014 data were all taken around phases 0.2–0.4 and show low values of p. On the other hand, the 2012 data were taken close to phase 0 and p is significantly higher, albeit with large uncertainties.

Zoom In Zoom Out Reset image size

We note, however, that the phase of the 2003 data is relatively uncertain, as the period was extrapolated back by over 10 years. With an estimated uncertainty of ±0.5 days on the period (W16), the 2003 data can be located at phases between −0.18 and 0.58. For the following discussion, we will nonetheless assume the predicted value of ∼0.2.

To check if the spectral variation with the 78-day phase is stable over longer periods of time, we extracted spectra from the Swift/XRT monitoring campaign between 2014 March and 2016 July using the online Swift/XRT data products generator (Evans et al. 2009). We selected data for the low-phase spectra between phases $0.7\leqslant \phi \leqslant 0.15$ and for the high phase between $0.3\leqslant \phi \leqslant 0.65$, based on a period of P = 78.12 days and ${T}_{0}=56663.0$ MJD, the minimum of the profile. While the data quality of these XRT spectra are not sufficient to confirm the spectral changes seen in the XMM-Newton data, they are fully consistent with the respective XMM-Newton models.

To look for possible degeneracies between p and temperature, we calculate confidence contours for these two parameters for each observation, shown in Figure 5. We find that the 2003 and 2014 data are compatible with each other, while both the 2012 and 2013 data are significantly different. The 2012 data are consistent within their uncertainties with a standard temperature profile for a thin disk (i.e., p = 0.75, Shakura & Sunyaev 1973).

Zoom In Zoom Out Reset image size

The uncertainties in p are also influenced by the variability of the absorption column. If we fix the column at $6.7\times {10}^{21}$ cm⁻², the average value of all four epochs, we find drastically smaller confidence contours while finding statistically comparable fits (${\rm{\Delta }}{\chi }^{2}\leqslant 3$ in all cases). All contours are fully within the contours presented in Figure 5; however, we note that the 2012 data are no longer consistent with values of $p\geqslant 0.75$, thus ruling out a standard thin disk. Because we do not know if the intrinsic absorption column has changed between observations, we continue our discussion with the model in which the column is allowed to vary independently between observations.

We also show contours for 2013 and 2014 using the XMM-Newton data only, to investigate if only using data below 10 keV results in a systematic parameter shift. These contours are shown in light color in Figure 5. We find that the XMM-Newton-only data seem to prefer a higher temperature, but find almost the same values for p. This is not surprising as the peak of the spectrum, which determines the temperature, is at the edge of the XMM-Newton range, while p is most relevant at lower energies and determines the slope of the spectrum up to the peak. Therefore, the XMM-Newton data drive p, while the temperature is mainly determined through the high-energy coverage of NuSTAR.

It is currently not known how the accretion disk is structured in neutron star ULXs. For the often implied strong magnetic fields of these systems (e.g., Dall'Osso et al. 2015; Mushtukov et al. 2015; Tong 2015), a truncation of the disk at ∼1400 km radii is expected. However, due to the very high luminosity of NGC 5907 ULX1, the spherization radius, where the thickness of the disk becomes comparable to the distance for the compact object (Shakura & Sunyaev 1973), is still about an order of magnitude larger than the magnetospheric radius (King & Lasota 2016). This implies that a geometrically thick disk is present, in which advection becomes important, in agreement with our spectral models. However, to sustain super-Eddington accretion rates, Israel et al. (2016a) argue that the disk has to be thin at the magnetospheric radius to not engulf the neutron star up to high magnetic latitudes. This problem could be mitigated by strong geometrical collimation, allowing the radiation to escape along a narrow funnel, which, however, is at odds with the observed smooth sinusoidal pulse profile (see also Bachetti et al. 2014; Fürst et al. 2016b).

We put the 2014 data into context with observations taken in 2003, 2012, and 2013. We find that in all observations the modeled disk temperature is similar, with ${T}_{\mathrm{in}}\approx 2.75\,\mathrm{keV}$. Assuming a slim-disk accretion geometry, we might therefore assume that we always see the hottest, innermost part of the accretion disk, which is likely located close to the magnetospheric radius ${r}_{{\rm{m}}}$. Israel et al. (2016a) argue that the surface magnetic field of NGC 5907 ULX1 should be around 3 × 10¹³ G, which implies ${r}_{{\rm{m}}}\approx 1400$ km (Cui 1997). If we assume that the observed radiation is dominated by the accretion disk, we can estimate the viewing angle from the normalization of the diskpbb model.

We follow the description of Soria et al. (2015) and Brightman et al. (2016). We assume a color correction factor of $\kappa =3$ and a geometric factor of $\xi =0.353$, appropriate for the high Eddington fraction of NGC 5907 ULX1 (Soria et al. 2015). We then calculate the inclination θ as

$$\begin{eqnarray}&&\cos (\theta )=\displaystyle \frac{{\xi }^{2}{\kappa }^{4}N{{d}_{10}}^{2}}{{{r}_{{\rm{m}}}}^{2}}.\end{eqnarray} \tag{ 1 }$$

Where N is the normalization of the diskpbb model and d₁₀ is the distance in units of 10 kpc. This results in an almost edge-on view of $\sim 89^\circ$.

This inclination is very high and makes it difficult to explain how we would be able to see the regions close to the compact object at all, as in super-Eddington accretion the accretion disk is expected to have some geometrical thickness, blocking our line of sight. The result depends on our assumptions of κ and ξ, but we obtain high inclinations for all realistic assumptions (e.g., $1.7\leqslant \kappa \leqslant 3$, Watarai & Mineshige 2003). About 20% of the observed flux is pulsed (Israel et al. 2016a), indicating that it is produced by the rotating accretion column of the neutron star. A lower flux from the disk, however, increases the estimated inclination further.

Another large systematic uncertainty in this estimate is the strength of the magnetic field, and the location of the hottest part of the disk. King & Lasota (2016) argue that a low magnetic field (∼10¹⁰ G) can explain the observed timing properties of M82 X-2, which are very similar to NGC 5907 ULX1. With a lower magnetic field and a consequently much smaller magnetospheric radius, the implied inclination would be smaller. In fact, in a face-on geometry ($\cos (\theta )=1$), we would imply a magnetic field of $B\approx 6\times {10}^{10}$ G for the inner disk radius to be at the magnetospheric radius during the bright 2014 observation.

The 2013 observation shows a distinctly different spectrum than the 2012 observation, despite having a very similar flux. It also shows a significantly different spectrum than both the 2003 and 2014 observations, despite being located in phase between them. This behavior might be related to the fact that the 2013 data were obtained only four days after the source's luminosity was below the detection limit of XMM-Newton and NuSTAR (W15). This remarkable drop in luminosity is clearly not related to the stable 78-day period, and so far has not been seen to repeat (W16). It is reasonable to expect that it is caused by the so-called "propeller effect" or the centrifugal inhibition of accretion, where the magnetospheric radius becomes larger than the corotation radius. This regime can be entered if the ram pressure of the accreting material drops, leading to a further dramatic reduction in accretion rate and consequently luminosity. During such a state, the inner accretion disk would get depleted, and the slightly lower temperature measured could be an indication that the inner accretion disk was still in the process of refilling during the second 2013 observation in which the source was detected.

W16 argued that the 78-day period is most likely either orbital or super-orbital in nature. The new results by Israel et al. (2016a) indicate an orbital period of ∼5 days, confirming the super-orbital nature of the 78-day period.

The physical origin of super-orbital periods is often linked to precession of the accretion disk (e.g., in Her X-1 and SMC X-1; Schandl & Meyer 1994; Caproni et al. 2006), though other possibilities have also been discussed (for a review, see Kotze & Charles 2012).

A precessing disk could explain the observed regular flux variations without changes in the physical conditions of the accretion flow because our viewing angle would change periodically (as also postulated for ULX M82 X-2, Kong et al. 2016). This model would not apply to the 2013 data because we very likely observed the source during an unusual state, as discussed above. Using the ratio of the diskpbb normalizations and the geometric effect that the flux is reduced by $\sqrt{\cos (\theta )}$, under the assumption that the disk is relatively flat, we can calculate the required change in inclination angle θ. The largest variance in θ is required if we observe NGC 5907 ULX1 face-on ($\theta =0$) during the bright phase. This would require an inclination angle of $\theta \approx 25^\circ$ during the faint phases (based on the 2012 data) and a half opening angle of the precession of $\sim 13^\circ$. This is lower than the precession seen in SS 433 (Khabibullin & Sazonov 2016), which is often argued to be a Galactic example of super-Eddington accretion analogous to ULXs, but viewed close to edge-on such that the X-ray emitting regions are obscured from view (Fabrika 2004). Under the assumption of a strong magnetic field, we estimated much higher viewing angles (almost edge-on, see Section 4.1), which would require a much smaller variability in θ to obtain similar flux variability. A variable viewing angle can therefore explain the observed flux changes for all inclinations (with the exception of the unusual 2013 data).

While a precessing disk can naturally explain the differences in observed flux across the 78-day cycle, we need to understand how different spectral shapes can be measured at different viewing angles. From Figure 4 it is clear that the strongest spectral change is observed in p. At super-Eddington accretion rates, it is expected that the accretion disk is flared up due to radiation pressure, i.e., it becomes geometrically thick and increases in size with radius. In addition, a strong wind is launched, which is also largely optically thick (Poutanen et al. 2007; Dotan & Shaviv 2011) and for which observational evidence has recently been found in NGC 1313 X-1 and NGC 5408 X-1 (Pinto et al. 2016; Walton et al. 2016b). The observed temperature profile, therefore, depends on which parts and with what angle we observe the accretion disk. Qualitatively, we can envision a geometry where the apparent temperature gradient in the disk is changing as a function of viewing angle. Detailed calculations of this model are, however, beyond the scope of this paper.

Israel et al. (2016a) did not find pulsations in the 2012 and 2013 XMM-Newton observations and give an upper limit of a pulsed fraction of 12%. Both of these observations were taken at low apparent luminosities. If these changes are connected to real changes in accretion rate, the properties of the accretion column might change. For a lower accretion rate, as observed in 2013, the emission pattern of the accretion column might be wider, resulting in a reduced pulsed fraction.

We have argued that the low observed flux of the 2012 data is not due to a lower intrinsic flux, but due to a change in viewing angle. This seems at first difficult to reconcile with the disappearing of pulsations. However, it is also possible that the neutron star shows free precession, in step with the precession of the accretion disk (e.g., as discussed for Her X-1, Staubert et al. 2013). In this case, the rotational axis at early super-orbital phases might be aligned close to our line of sight or, if the emission is axial-symmetric, close to 90°, also reducing the observable pulsed fraction.

It is clear that, to understand super-Eddington neutron stars like NGC 5907 ULX1, more observational and theoretical work needs to be done. For example, the currently available coverage of the phase-space is concentrated between phases 0–0.5. For a detailed test of the proposed spectral evolution and model, broadband observations at later phases are necessary, which can be obtained with XMM-Newton and NuSTAR. NGC 5907 ULX1 is an ideal target for these studies, given the stable and strong 78-day period. Its spectral similarities to other ULXs will help us understand this class of objects better and investigate how the type of compact object influences their behavior.