Super-Eddington Accretion onto the Neutron Star NGC 7793 P13: Broadband X-ray Spectroscopy and Ultraluminous X-ray Sources

TL;DR

This study provides a detailed broadband X-ray spectral analysis of the ULX pulsar NGC 7793 P13, revealing complex accretion flow structures and magnetic field constraints, and compares its spectral features with other ULXs to suggest neutron star accretors may be common.

Contribution

It offers the first broadband spectral analysis of P13, identifies multiple thermal components, and constrains the neutron star's magnetic field, advancing understanding of super-Eddington accretion in ULXs.

Findings

Two thermal blackbody components with distinct temperatures are required to fit the data.

The accretion flow outside the magnetosphere is complex, indicating a possible truncated inner disc.

The neutron star's magnetic field is constrained to be less than approximately 6 x 10^{12} G.

Abstract

We present a detailed, broadband X-ray spectral analysis of the ULX pulsar NGC 7793 P13, a known super-Eddington source, utilizing data from the $XMM$-$Newton$, $NuSTAR$ and $Chandra$ observatories. The broadband $XMM$-$Newton+NuSTAR$ spectrum of P13 is qualitatively similar to the rest of the ULX sample with broadband coverage, suggesting that additional ULXs in the known population may host neutron star accretors. Through time-averaged, phase-resolved and multi-epoch studies, we find that two non-pulsed thermal blackbody components with temperatures $\sim$0.5 and $\sim$1.5 keV are required to fit the data below 10 keV, in addition to a third continuum component which extends to higher energies and is associated with the pulsed emission from the accretion column. The characteristic radii of the thermal components appear to be similar, and are too large to be associated with the neutron star itself, so the need for two components likely indicates the accretion flow outside the magnetosphere is complex. We suggest a scenario in which the thick inner disc expected for super-Eddington accretion begins to form, but is terminated by the neutron star's magnetic field soon after its onset, implying a limit of $B \lesssim 6 \times 10^{12}$ G for the dipolar component of the central neutron star's magnetic field. Evidence of similar termination of the disc in other sources may offer a further means of identifying additional neutron star ULXs. Finally, we examine the spectrum exhibited by P13 during one of its unusual 'off' states. These data require both a hard powerlaw component, suggesting residual accretion onto the neutron star, and emission from a thermal plasma, which we argue is likely associated with the P13 system.

Tables

Table 1. Table 1: Details of the X-ray observations of NGC 7793 P13 considered in this work

Epoch	Mission(s)	OBSID(s)	Start	Exposure\tmark[a]
			Date	(ks)
2016
XN1	NuSTAR	80201010002	2016-05-20	115
	XMM-Newton	0781800101		28/46
Archival Data
C1	Chandra	3954	2003-09-06	49
C2	Chandra	14231	2011-08-13	59
C3	Chandra	13439	2011-12-25	58
C4	Chandra	14378	2011-12-30	25
X1	XMM-Newton	0693760101	2012-05-14	26/36
X2	XMM-Newton	0693760401	2013-11-25	41/47
X3	XMM-Newton	0748390901	2014-12-10	42/48

Table 2. Table 2: Best fit parameters obtained for the broadband continuum models applied to the P13 data from epoch XN1

Model\tmark[a]	$N_{H;\mathrm{int}}$	$k{T}_{\mathrm{in};\mathrm{DBB}}$	Norm	$p/\tau /{\Gamma }_{\mathrm{CPL}}$	$kT\text{tmark}\left[b\right]/E_{\mathrm{cut}}$	Norm	${\Gamma }_{\mathrm{SIMPL}}$	$f_{\mathrm{scat}}$	${\chi }^2$/DoF
	[${10}^{20}$ cm-2]	[keV]			[keV]			[%]
DISKPBB	$8\pm 1$	${0.45}_{-0.04}^{+0.03}$	${1.1}_{-0.2}^{+0.4}$	$>1.3$\tmark[c]	$2.0\pm 0.3$	$8_{-5}^{+6}\times {10}^{-2}$	${3.7}_{-1.1}^{+0.9}$	$>19$	1129/1148
COMPTT	$8\pm 1$	${0.49}_{-0.20}^{+0.10}$	${1.0}_{-0.5}^{+2.3}$	$17\pm 2$	${1.4}_{-0.1}^{+0.3}$	${7.0}_{-1.0}^{+0.6}\times {10}^{-4}$	${4.0}_{-0.7}^{+0.3}$	$>36$	1129/1148
CUTOFFPL	$8\pm 1$	${0.48}_{-0.05}^{+0.08}$	$1.0\pm 0.4$	$-{0.9}_{-0.6}^{+0.4}$	${1.9}_{-0.6}^{+0.7}$	$\left(2.5\pm 0.5\right)\times {10}^{-4}$	${4.0}_{-1.8}^{+0.7}$	$>12$	1129/1148

Table 3. Table 3: Best fit parameters obtained from our phase-resolved analysis of epoch XN1

Component	Parameter
TBABS	$N_{H;\mathrm{int}}$	[${10}^{20}$ cm-2]	$8_{-1}^{+2}$
DISKBB	$k{T}_{\mathrm{in}}$	[keV]	${0.40}_{-0.07}^{+0.01}$
	Norm		${1.5}_{-0.1}^{+1.4}$
DISKPBB	$k{T}_{\mathrm{in}}$	[keV]	$1.9\pm 0.2$
	$p$		$>0.95$\tmark[a]
	Norm	[${10}^{-2}$]	${4.6}_{-3.6}^{+0.8}$
CUTOFFPL	$\Gamma$		$-0.1$\tmark[b]
	$E_{\mathrm{cut}}$	[keV]	$3.9$\tmark[b]
	$F_{2-10}$\tmark[c] (low)	[${10}^{-12}$ erg cm-2 s-1]	${1.7}_{-0.3}^{+0.2}$
	$F_{2-10}$\tmark[c] (med)	[${10}^{-12}$ erg cm-2 s-1]	${2.6}_{-0.3}^{+0.2}$
	$F_{2-10}$\tmark[c] (high)	[${10}^{-12}$ erg cm-2 s-1]	${3.6}_{-0.3}^{+0.2}$
${\chi }^2$/DoF			942/932

Table 4. Table 4: Evolution in the Properties of the X-ray Pulsations across the Various Epochs Considered

Epoch	$F_{\mathrm{av},0.3-10}$	$P$	$P_{\mathrm{frac},0.5-2.0}$	$P_{\mathrm{frac},2-10}$	$F_{\mathrm{pulse},2-10}$
	[${10}^{-12}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$]	[ms]	[%]	[%]	[${10}^{-12}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$]
X2	$1.14\pm 0.03$	$419.712\pm 0.008$	$12\pm 2$	$31\pm 3$	$0.57\pm 0.08$
X3	$2.84\pm 0.05$	$418.390\pm 0.008$	$18\pm 2$	$29\pm 2$	$1.19\pm 0.12$
XN1	$5.19\pm 0.07$	$416.9513\pm 0.0017$	$13\pm 1$	$22\pm 1$	$1.84\pm 0.12$

Table 5. Table 5: Results obtained for the multi-epoch analysis of P13

Model	Parameter		Global	Epoch
Component				X2	X3	XN1
Continuum: DISKBB+DISKPBB+CUTOFFPL
TBABS	$N_{H,\mathrm{int}}$	[${10}^{20}$ cm-2]	$7.7\pm 0.8$
DISKBB	$T_{\mathrm{in}}$	[keV]		$0.34\pm 0.01$	${0.417}_{-0.026}^{+0.008}$	${0.421}_{-0.022}^{+0.009}$
	Norm			${1.24}_{-0.10}^{+0.32}$	${0.85}_{-0.07}^{+0.25}$	${1.25}_{-0.11}^{+0.25}$
DISKPBB	$T_{\mathrm{in}}$	[keV]		$1.70\pm 0.15$	${2.14}_{-0.11}^{+0.14}$	${2.15}_{-0.06}^{+0.10}$
	$p$		$>1.52$\tmark[a]
	Norm	[${10}^{-2}$]		${1.3}_{-0.5}^{+0.4}$	${1.9}_{-0.7}^{+0.3}$	${4.0}_{-1.4}^{+0.1}$
CUTOFFPL	$\Gamma$		$0.17$\tmark[b]
	$E_{\mathrm{cut}}$	[keV]	$4.7$\tmark[b]
	$F_{2-10}$	[${10}^{-13}$ erg cm-2 s-1]		$5.5$\tmark[c]	$11.3$\tmark[c]	${17.6}_{-2.8}^{+2.0}$
${\chi }^2$/DoF			2149/2058
Continuum: DISKBB+BB+CUTOFFPL
TBABS	$N_{H,\mathrm{int}}$	[${10}^{20}$ cm-2]	${7.1}_{-0.7}^{+0.8}$
DISKBB	$T_{\mathrm{in}}$	[keV]		$0.35\pm 0.02$	$0.45\pm 0.03$	$0.48\pm 0.03$
	Norm			${1.05}_{-0.21}^{+0.27}$	${0.63}_{-0.13}^{+0.17}$	${0.80}_{-0.15}^{+0.19}$
BB	$T_{\mathrm{in}}$	[keV]		$1.1\pm 0.1$	$1.5\pm 0.1$	$1.53\pm 0.04$
	Norm	[${10}^{-6}$]		$2.4\pm 0.6$	$11.1\pm 1.5$	${27.3}_{-2.2}^{+2.4}$
CUTOFFPL	$\Gamma$		$0.17$\tmark[b]
	$E_{\mathrm{cut}}$	[keV]	$4.7$\tmark[b]
	$F_{2-10}$	[${10}^{-13}$ erg cm-2 s-1]		$6.9$\tmark[c]	$14.4$\tmark[c]	${22.3}_{-1.6}^{+1.4}$
${\chi }^2$/DoF			2144/2059

Table 6. Table 6: Best fit parameters obtained for the off-state datasets

Parameter		Global	Chandra	XMM
$k{T}_{\mathrm{mekal}}$	[keV]	${1.3}_{-0.2}^{+0.3}$
Norm	[${10}^{-6}$]	$5\pm 2$
$\Gamma$		${0.7}_{-0.5}^{+0.4}$
Norm	[${10}^{-6}$]		${1.8}_{-0.9}^{+1.2}$	${0.8}_{-0.5}^{+0.9}$
${\chi }^2$/DoF		28/27
$L_{\mathrm{tot}}$\tmark[a]	[${10}^{37}$ erg s-1]		${7.9}_{-1.3}^{+1.2}$	${4.5}_{-1.6}^{+1.4}$
$L_{\mathrm{po}}$\tmark[a]			$6.5\pm 1.4$	${3.0}_{-1.7}^{+1.2}$

Full text

††pagerange: Super-Eddington Accretion onto the Neutron Star NGC 7793 P13: Broadband X-ray Spectroscopy and Ultraluminous X-ray Sources–References

1 Introduction

The ultraluminous X-ray source (ULX111Here we define a ULX to be an off-nuclear X-ray point source with a luminosity in excess of $10^{39}$ erg s*-1* in the 0.3–10.0 keV bandpass, assuming isotropy.) NGC 7793 P13 (hereafter P13), which reaches luminosities of $L_{\rm{X,peak}}\sim 10^{40}$ erg s*-1*, was recently found to be powered by an accreting neutron star through the detection of coherent X-ray pulsations (Fürst et al. 2016b; Israel et al. 2017b). Along with M82 X-2, the first ULX pulsar discovered ($L_{\rm{X,peak}}\sim 2\times 10^{40}$ erg s*-1*; Bachetti et al. 2014) and NGC 5907 ULX ($L_{\rm{X,peak}}\sim 9\times 10^{40}$ erg s*-1*; Israel et al. 2017a; Fürst et al. 2017), only three neutron star ULXs with such extreme luminosities are currently known. These systems are of great interest, as they offer a rare opportunity to study accretion in the highly super-Eddington regime (the Eddington limit for a standard 1.4 $M_{\odot}$ neutron star is $L_{\rm{E}}\sim 2\times 10^{38}$ erg s*-1*).

There are several observational similarities seen between the P13, M82 X-2 and NGC 5907 ULX systems. All three have broadly similar pulse periods, of order $\sim$1 s (Bachetti et al. 2014; Fürst et al. 2016b; Israel et al. 2017a, b), and exhibit strong secular spin-up owing to the large accretion torques related to their extreme accretion rates; over just $\sim$10 years, the pulse period in NGC 5907 ULX evolved from $\sim$1.4 to $\sim$1.1 s. The pulse profiles observed are also similar, with all three sources showing broad, near-sinusoidal variations. Furthermore, all three show long-term periodicities ($\sim$60–80 d; Motch et al. 2014; Walton et al. 2016a; Kong et al. 2016; Hu et al. 2017, Brightman et al. in prep.). In the cases of M82 X-2 and NGC 5907 ULX these periods are known to be super-orbital, as variations in their pulse periods have also revealed orbital periods of a few days. No such independent constraints on the orbit are currently available for P13, where the nature of the $\sim$64 d period remains uncertain, although it does seem to show some ‘jitter’ (Motch et al. 2014) that may be indicative of super-orbital variations. Finally, all three sources show unusual ‘off-states’, where the X-ray flux is diminished by a factor of $\sim$50 or more (e.g. Motch et al. 2014; Walton et al. 2015a; Brightman et al. 2016a). Tsygankov et al. (2016) suggest that these extreme flux modulations may be related to the propeller effect, but their nature is not currently well established, and may even differ between systems.

These systems pose key questions, particularly with regards to how these neutron stars are able to reach such extreme apparent luminosities. The magnetic field of the central neutron star will channel the accretion flow into columns inside the magnetospheric radius ($R_{\rm{M}}$, the point at which magnetic pressure dominates), allowing material to accrete onto the magnetic poles while radiation escapes from the sides of the column (Basko & Sunyaev, 1976). This introduces anisotropy into the radiation field (which is required in order for X-ray pulsations to be observed), and in turn a beaming correction to simple luminosity estimates becomes necessary. However, while some anisotropy is clearly required, it is difficult to explain the near-sinusoidal pulse profiles observed from all three systems in the context of strong beaming of sub-Eddington accretion onto an otherwise standard neutron star (which would require beaming by factors of $\sim$100 or more).

Some authors have invoked strong, magnetar-level magnetic fields ($B\sim 10^{13-14}$ G; e.g. Ekşi et al. 2015; Dall’Osso et al. 2015; Mushtukov et al. 2015). This reduces the scattering cross section for electrons (Herold, 1979), reducing the radiation pressure and in turn increasing the effective Eddington luminosity. However, other authors have instead suggested that the magnetic field is much lower (potentially as low as $B\sim 10^{9}$ G) based on the ratio of the spin-up rate to the luminosity, which is an order of magnitude lower than typical X-ray pulsars (e.g. Kluźniak & Lasota 2015). These authors argue that a disk truncated at a large radius (close to the co-rotation radius, $R_{\rm{co}}$, as would occur for a high-B field system; $R_{\rm{co}}$ is the point at which the material in the disc co-rotates with the neutron star, and $R_{\rm{M}}$ $<$ $R_{\rm{co}}$ is required for strong accretion to occur) would not provide the required lever arm to power the observed spin up. In such a scenario, the extreme luminosities would need to be produced by a highly super-Eddington accretion disk that extends close to the accretor, similar to super-Eddington accretion onto a black hole (King & Lasota 2016).

P13 provides a key laboratory for understanding these extreme neutron stars and their link to the broader ULX population. With only three examples, the sample of ULX pulsars is still extremely limited, and there are obstacles to undertaking detailed analyses with M82 X-2 and NGC 5907 ULX. M82 X-2 is heavily confused with its more luminous neighbour M82 X-1 (separated by $\sim$5*′′), and NGC 5907 ULX is situated at a large distance ($\sim$17 Mpc; Tully et al. 2016) so its observed flux is relatively low. Both sources are also heavily absorbed ($N_{\rm{H}}\sim 10^{22}$ cm-2*), severely complicating comparisons with other ULXs. In contrast, P13 is both well isolated and relatively unobscured ($N_{\rm{H}}\sim 10^{21}$ cm*-2*). Furthermore, pulsations are only intermittently detected in M82 X-2 and NGC 5907 ULX (Bachetti et al. 2014; Israel et al. 2017a), while for P13 they have always been detected when the source has been observed at high flux and with sufficient temporal resolution.

In this paper, we present results from the 2016 broadband observation of P13 with XMM-Newton (Jansen et al. 2001) and NuSTAR (Harrison et al. 2013). Our initial work on these data, presented in Fürst et al. (2016b), focused on the detection and evolution of the X-ray pulsations, and here we focus on undertaking a comprehensive spectral analysis. We describe the observations and data reduction in Section 2, and we provide details on our spectral analysis in Sections 3 and 4. We draw comparisons with the broader ULX population in Section 5, and in Section 6 we discuss the results. Finally, we summarise our conclusions in Section 7. Throughout this work, we adopt a distance to NGC 7793 of 3.5 Mpc (Pietrzyński et al. 2010).

2 Observations and Data Reduction

As discussed in Fürst et al. (2016b), the coordinated XMM-Newton+NuSTAR observation of P13 that is the focus of much of this work was triggered following the detection of a high flux with Swift (Gehrels et al. 2004), and currently represents the highest flux probed with high signal-to-noise (S/N) observations for this source to date. In addition to these data, we also consider the archival XMM-Newton and Chandra (Weisskopf et al. 2002) exposures of P13. Details of all the observations considered in this work are given in Table 1. Our data reduction procedure largely follows that outlined in Fürst et al. (2016b), but here we repeat the salient points for clarity.

2.1 NuSTAR

The NuSTAR data were reduced using the NuSTAR Data Analysis Software (NUSTARDAS; v1.6.0) following standard procedures; NuSTAR caldb v20160824 was used throughout. We cleaned the unfiltered event files with NUPIPELINE, applying the standard depth correction (which significantly reduces the internal background) and removing passages through the South Atlantic Anomaly. Barycentre corrections were applied to the photon arrival times in the cleaned event files using using the DE200 reference frame, and based on the pulse period reported in Fürst et al. (2016b) we phase-stamp each of the events. Source products and instrumental response files were extracted from circular regions of radius 70*′′* for both of the focal plane modules (FPMA/B) with NUPRODUCTS, with background estimated from a much larger region on the same detector as P13. In addition to the standard ‘science’ data, to maximize the signal-to-noise (S/N) we also extract the ‘spacecraft science’ data following the procedure outlined in Walton et al. (2016c), which provides $\sim$10% of the good exposure in this case.

2.2 XMM-Newton

The XMM-Newton data were reduced with the XMM-Newton Science Analysis System (v15.0.0), following the standard data reduction guide.222http://xmm.esac.esa.int/. The raw data files were cleaned using EPCHAIN and EMCHAIN for the EPIC-pn (Strüder et al. 2001) and each of the EPIC-MOS (Turner et al. 2001) detectors, respectively. As with the NuSTAR data, barycentre corrections were applied to the photon arrival times in the cleaned event files for EPIC-pn using the DE200 reference frame, and we phase-stamp events based on the pulse periods reported in Fürst et al. (2016b); unfortunately the EPIC-MOS detectors do not have sufficient temporal capabilities to resolve the $\sim$0.4 s pulse period. Source products were typically extracted from circular regions of radius $\sim$40*′′* and background was estimated from larger areas on the same CCD free of other contaminating point sources. The exception was OBSID 0693760101, during which the source was extremely faint (but still clearly detected), so a radius of 20*′′* was used. As recommended, we only use single–double patterned events for EPIC-pn and single–quadruple patterned events for EPIC-MOS, and we also excluded periods of high background as standard. Instrumental redistribution matrices and ancillary response files were generated for each of the detectors with RMFGEN and ARFGEN, respectively. Finally, after checking their consistency, spectra from the two EPIC-MOS detectors were combined using ADDASCASPEC.

2.3 Chandra

The Chandra data were reduced with the CIAO software package (v4.7) following standard procedure. We extracted spectra from the ACIS-S detector (Garmire et al. 2003) level 2 event files and computed instrumental response files using SPECEXTRACT. Although the three 2011 observations all occurred during the extended off-state seen from P13, the source is still clearly detected in all observations. Source spectra are extracted from circular regions of radius 4–5*′′*, and as with both NuSTAR and XMM-Newton background was estimated from larger areas on the same CCD free of other contaminating point sources. All the observations were performed in the Timed Exposure mode with fixed frame times of $\sim$3.2 s; the time-resolution of these Chandra observations is therefore too coarse to detect the $\sim$0.4 s pulsations from P13.

3 NGC 7793 P13: The Broadband 2016 Dataset

We begin our analysis by focusing on the broadband XMM-Newton+NuSTAR dataset obtained in 2016 (*i.e. *epoch XN1), as these observations have provided the first ever hard X-ray ($>$10 keV) detection of P13. Model fits are performed with XSPEC  v12.6.0f (Arnaud 1996), parameter uncertainties are quoted at the 90% confidence level for one interesting parameter ($\Delta\chi^{2}=2.7$, with all other free parameters allowed to vary during the error estimation). All models considered include a Galactic absorption component with a fixed column of $N_{\rm{H,Gal}}=1.2\times 10^{20}$ $\rm\thinspace cm^{-2}$ (Kalberla et al. 2005), and we also allow for absorption intrinsic to the source at the redshift of NGC 7793 ($N_{\rm{H;int}}$; $z=0.000767$). Both neutral absorption components are modelled with the TBNEW absorption code. As recommended for TBNEW, we use the cross-sections of Verner et al. (1996) and the abundance set presented by Wilms et al. (2000). We also allow for cross-calibration uncertainties between the different detectors by including multiplicative constants that are allowed to float between the datasets, fixing EPIC-pn at unity. These constants are always within 10% of unity, as expected (Madsen et al. 2015). Unless stated otherwise, spectra are grouped to have a minimum of 50 counts per energy bin throughout this work in order to facilitate the use of $\chi^{2}$ minimization during our analysis.

3.1 Time-Averaged Spectroscopy

The time-averaged broadband spectrum observed from epoch XN1 is shown in Figure 1 (left panel); P13 is detected by NuSTAR up to $\sim$30 keV. The broadband spectrum is qualitatively similar to the broadband spectra seen from the rest of the ULX population to date (see Section 5). As noted by Motch et al. (2014) and Pintore et al. (2017), who present spectral analyses of the archival XMM-Newton observations of P13, the spectrum below 10 keV requires two broadband continuum components: a thermal component that contributes below $\sim$2 keV, and a harder component that dominates at higher energies. The prior XMM-Newton observations analyzed by these authors indicate this higher energy component shows curvature in the $\sim$5–10 keV band, similar to the rest of the ULX population (e.g. Stobbart et al. 2006; Gladstone et al. 2009). The NuSTAR data clearly and robustly confirms the presence of this curvature, demonstrating that, as with the rest of the ULX sample observed by NuSTAR to date, the emission that dominates the 2-10 keV band is not a high-energy powerlaw tail, and that the curvature seen in the XMM-Newton data does mark the start of a high-energy spectral cutoff.

Before undertaking a more detailed analysis of the spectral properties of P13, we fit the average spectrum with a set of models typically applied to other ULXs in the literature to provide a simple comparison. Based on the previous analyses of the archival P13 data, and our prior work on broadband observations of ULXs with NuSTAR, we apply a set of 2-component continuum models commonly applied to ULX spectra. We start with a model consisting of two accretion disk components, combining the DISKBB and DISKPBB models (Mitsuda et al. 1984; Mineshige et al. 1994). The former is based on the thin accretion disc model of Shakura & Sunyaev (1973), which has a radial temperature index of $p=3/4$ (where $T(r)\propto r^{-p}$) and just has the disc temperature ($T_{\rm{in}}$) as a free parameter, while the latter also allows $p$ to be varied as a free parameter. Second, we also apply a model combining DISKBB with COMPTT (Titarchuk 1994), a thermal Comptonization model which is primarily characterized by the optical depth ($\tau$) and the temperature ($kT_{\rm{e}}$) of the scattering electrons (for convenience, we assume the temperature of the seed photons to be that of the DISKBB component, as in previous works). Both of these models provide reasonable fits to the data below $\sim$10–15 keV. However, as can clearly be seen in the data/model ratios shown in Figure 1 (right panels), both models leave a clear excess in the NuSTAR data at the highest energies probed.

To further test the presence of this hard excess, we also fit a third model, combining DISKBB with a more phenomenological CUTOFFPL component. This is a simple powerlaw continuum with an exponential cutoff ($E_{\rm{cut}}$), which has a broader curvature at high energies than both of the DISKPBB and COMPTT models discussed above. This model has often been applied to ULX data (e.g.Bachetti et al. 2013; Pintore et al. 2017), but is also fairly representative of models typically applied to Galactic neutron stars (e.g. Coburn et al. 2002; Fürst et al. 2014)333We have also tried replacing the CUTOFFPL model with the FDCUT and NPEX models also frequently applied to Galactic neutron stars (e.g. Odaka et al. 2013; Fürst et al. 2013; Islam et al. 2015; Vybornov et al. 2017), but these provide identical fits to the CUTOFFPL model, so we do not present the results in full.. However, although it is weaker than both the previous models, even the broader curvature cannot resolve the hard excess (see Figure 1); a third continuum component is required to fit the time-averaged data regardless of the model used to fit the curvature in the $\sim$5–10 keV band. Similar hard excesses have now been seen in a number of other ULXs observed with NuSTAR (e.g. Walton et al. 2015b, 2017; Mukherjee et al. 2015; Fürst et al. 2017). Following these works, to account for this excess emission we add SIMPL (Steiner et al. 2009) to the higher energy component in the models described above. This is a convolution model that scatters some fraction ($f_{\rm{scat}}$) of the photons in an input continuum model into a high-energy powerlaw tail (*i.e. *photon number is conserved). The additional of SIMPL resolves the hard excess and results in a significant improvement to the fit for all three of the base continuum models considered; even for the DISKBB+CUTOFFPL model – for which the hard excess is weakest – the fit improves by $\Delta\chi^{2}=32$ for two additional free parameters.

In order to confirm that this additional high-energy component is required by the data at a significant level we performed a series of simulations, using the DISKBB+CUTOFFPL continuum model to be conservative, as this is the model that returns the lowest $\Delta\chi^{2}$ when this component is included in the fit. Using the same response and background files, and adopting the same exposure times as the real data, we simulated 10,000 sets of XMM-Newton (pn and combined MOS1+MOS2) and NuSTAR (FPMA and FPMB) spectra with the FAKEIT command in XSPEC based on the best-fit DISKBB+CUTOFFPL model. Each of the simulated datasets was rebinned in the same manner and analysed over the same bandpass as the real data. We then fit each of the simulated XMM-Newton+NuSTAR datasets with the DISKBB+CUTOFFPL, and noted the $\Delta\chi^{2}$ improvement the addition of a SIMPL component provided over this fit. Of the 10,000 datasets simulated, none returned a chance improvement equivalent to or greater than that observed, implying that this component is seen in the real data with a significance comfortably in excess of 3$\sigma$ (for reference, the false-alarm probability for a 3$\sigma$ detection would correspond to $\Delta\chi^{2}\sim 9$ here).

The results for the three models considered here, which all provide statistically equivalent fits, are presented in Table 2. While the broadband spectrum is qualitatively similar to the rest of the sample of ULXs observed by NuSTAR to date, there are some notable quantitative differences. For example, the temperature of the lower energy DISKBB component is $kT_{\rm{in;DBB}}\sim 0.5$ keV, slightly hotter than for other ULXs where the same models usually find $kT_{\rm{in;DBB}}\sim 0.3$ keV. Furthermore, the rise of the spectrum up to the peak of the emission at $\sim$7 keV is harder than typically observed from other ULXs. This can be seen in the results for the radial temperature index in the DISKPBB model and the optical depth in the COMPTT model; in the former case we find $p>1.3$ while ULXs more typically show $p<0.75$ when fit with DISKPBB, and in the latter case we find $\tau\sim 17$ while typical ULX values are $\tau\sim 5-10$ with this model. This means the continuum component dominating the 2-10 keV emission is much more peaked (*i.e. *the spectral curvature in this band is stronger) in P13 than in the broader ULX population. Indeed, if we replace this component with a single-temperature blackbody, a similarly good fit is obtained ($\chi^{2}$/DoF = 1133/1149), so a range of temperatures is not strongly required in this case.

3.2 The Pulsed Emission

To further characterize the spectral behaviour exhibited by P13, we also investigate spectral variations across the pulse period in an attempt to separate the emission from the accretion column (pulsed) and rest of the accretion flow (which for simplicity is assumed to be steady over the pulse cycle). Owing to its superior time resolution, in this section we are only able to utilize data from the XMM-Newton EPIC-pn detector for the lower energy data; the pulse period is too short to be resolved by the timing capabilities of the EPIC-MOS detectors.

3.2.1 Difference Spectroscopy

In order to undertake a comparison between P13 and M82 X-2, we compute the spectrum of the pulsed component for the broadband XMM-Newton+NuSTAR dataset. This approach is necessary since the absolute spectrum of M82 X-2 cannot easily be disentangled from the contribution from M82 X-1. We follow the approach taken for M82 X-2 in Brightman et al. (2016a), and subtract the average spectrum obtained during the faintest quarter (i.e. $\Delta\phi_{\rm{pulse}}=0.25$) of the pulse cycle from that extracted during the brightest quarter (i.e. “pulse on”$-$“pulse-off”) in order to isolate the pulsed emission in a simple manner. The data are rebinned to have a minimum S/N per energy bin of 3, to allow the use of $\chi^{2}$ statistics, and we fit the data over the 0.3–20.0 keV energy range with a simple CUTOFFPL model. The intrinsic column density is poorly constrained in these fits, and so we fix it to $8\times 10^{20}$ $\rm\thinspace cm^{-2}$ following our analysis of the average spectrum, and we also fix the cross-calibration constants between EPIC-pn and FPMA/B to the best-fit values obtained in that analysis given the low S/N of the pulsed spectrum. This provides a good fit to the data, with $\chi^{2}=113$ for 91 DoF. However, for P13 we find the pulsed emission to have a much harder rise, and a lower cutoff energy: $\Gamma=-0.1\pm 0.3$, $E_{\rm{cut}}=3.9^{+0.8}_{-0.6}$ keV. Identical results are obtained if we minimise the Cash statistic (Cash 1979) instead. For comparison, we show the 2-D confidence contours for these parameters in Figure 2 for both P13 and M82 X-2. The pulsed flux from P13 during this epoch corresponds to an apparent 0.3–20.0 keV luminosity of $(4.0\pm 0.3)\times 10^{39}$ erg s-1 were one to assume isotropic emission (although clearly this is not the case).

3.2.2 Phase-Resolved Spectroscopy

In addition to the “pulse-on”$-$“pulse-off” spectroscopy, we also fit the spectra extracted across a series of phase bins. We continue using phase bins of $\Delta\phi_{\rm{pulse}}=0.25$ in size, and extract spectra from four periods across the pulse cycle: peak, fall, trough and rise. However, on inspection the fall and rise spectra were found to be practically identical, so we combined them to form a higher S/N medium-flux spectrum, resulting in three broadband spectra to be fit in our phase-resolved analysis: high-flux (peak), medium-flux (rise $+$ fall) and low-flux (trough). We show the high- and low-flux spectra in Figure 3 for comparison. It is clear that the majority of the difference between the two is seen at higher energies, which is not surprising since the pulse fraction is known to increase with energy (Fürst et al. 2016b; Israel et al. 2017b); the observed spectrum is hardest during the peaks of the pulse cycle.

In order to model these data, we fit all three of the pulse-resolved phase bins simultaneously and undertake a simple decomposition of these spectra into stable and variable (*i.e. *pulsed) components. Following the previous section, we treat the pulsed emission simply with a CUTOFFPL model. The ‘shape’ parameters for this component (i.e. $\Gamma$, $E_{\rm{cut}}$) are linked across all the phase bins and fixed to the best-fit results from the difference spectroscopy described above, but its normalization can vary between them. We do not require the normalization of the pulsed CUTOFFPL component to be zero for the low-flux data, as the emission from the accretion column can still contribute during the minimum of the pulse cycle. For the stable emission, we return to the 2-component continuum models discussed previously, and apply the DISKBB+DISKPBB combination. All the parameters for these components, including their normalisations, are linked across all the phase bins, as is the intrinsic neutral absorption column.

The global fit to the phase-resolved data with this approach is excellent ($\chi^{2}$/DoF = 942/932), and the results are presented in Table 3. The results for the constant DISKBB+DISKPBB continuum components are broadly similar to the results obtained with the model utilizing this combination in our time-averaged spectral analysis (see Table 2). This suggests that for the case of P13, the need for the SIMPL component in the time-averaged models is largely driven by the contribution of the pulsed emission. The relative contribution of the DISKBB, DISKPBB and CUTOFFPL components in our analysis is shown for the peaks and the troughs of the pulse cycle in Figure 4; the CUTOFFPL model dominates the flux at the highest energies probed by NuSTAR, exactly where the SIMPL component contributes in the time-averaged models. Indeed, the steady continuum emission does not require an additional high-energy powerlaw tail in these fits; adding a SIMPL component to the DISKPBB component only improves the fit by $\Delta\chi^{2}$ = 6 for two additional free parameters (although the situations are not identical, this would be below the threshold required to claim a significant detection based on the simulations presented above). We stress though that both the DISKBB and DISKPBB components are required to describe the phase-resolved data; removing either significantly degrades the global fit. Three continuum components are therefore required to fully describe the phase-resolved data from epoch XN1: one pulsed component and two steady components.

4 NGC 7793 P13: Multi-Epoch analysis

We next undertake a multi-epoch analysis for P13 additionally including the archival XMM-Newton datasets (see Table 1). A comparison of the average spectra from each of these observations with the broadband XMM-Newton+NuSTAR dataset is shown in Figure 5. The first of the archival XMM-Newton observations (epoch X1) was taken during the extended off-state that spanned mid 2011 to mid 2013 (Motch et al. 2014; Fürst et al. 2016b) and subsequently only a very low S/N spectrum is available. However, the latter two both caught P13 at ULX luminosities, with fluxes only a factor of $\sim$4 and $\sim$2 fainter than epoch XN1 during epochs X2 and X3, respectively (see Table 4). For the ULX-luminosity observations, the 0.3–10.0 keV spectrum of P13 is clearly harder when the source is brighter, and exhibits stronger variations at the higher energies of this bandpass. During the faintest of these observations (epoch X2) the spectrum is most visibly doubly-peaked, as the relative contribution of the soft thermal component is strongest during this epoch.

4.1 Evolution of the Pulsed Emission

As noted by Fürst et al. (2016b), during epoch X1 P13 was too faint for any pulsations to be detectable. In addition, none of the Chandra observations have sufficient temporal resolution to see the pulsations. However, pulsations were detected in all the subsequent epochs (X2, X3 and XN1). Table 4 summarizes the results obtained from the pulsation searches presented in Fürst et al. (2016b), showing a clear spin-up. In addition to these results, we also compute the pulse fractions for the 0.5–2.0 and 2.0–10.0 keV bands that are common to all three epochs where pulsations have been detected, based on the XMM-Newton EPIC-pn lightcurves. Here, we calculate the pulse fraction as (max-min)/(max+min), where max and min refer to the average count rates at the maximum and minimum phases of the pulse cycle. The results are shown in Figure 6.

As expected, the pulse fractions are systematically higher in the harder 2.0–10.0 keV band. At these energies, the pulse fraction shows evidence for a secular decrease with increasing observed flux. This suggests that there are two emission components contributing in the 2.0–10.0 keV band (pulsed and non-pulsed), and that the relative contribution from the pulsed component is lower at higher fluxes. In contrast, in the softer band the pulse fraction initially increases with flux, and then decreases again at even higher fluxes. The evolution seen at these energies between epochs X3 and XN1 (the two higher-flux observations) is likely related to the same effect producing the anti-correlation in the harder band, *i.e. *the relative contribution of the pulsed component decreasing relative to some non-pulsed component as the flux increases. The switch to a decreasing pulse fraction with decreasing flux is then likely related to the increased relative importance of the cooler thermal component during epoch X2, which is visually apparent from Figure 5. This implies that the cooler thermal component is also not pulsed, resulting in the lower soft-band pulse fractions.

The relative evolutions seen in the pulse fractions for these two bands therefore suggests that there are two emission components that contribute to both the hard and the soft bands (pulsed and non-pulsed), and a third non-pulsed component that only contributes to the soft band. This is broadly similar to the spectral decomposition presented in Section 3.2.

We also investigate the spectral evolution of the pulsed emission, by repeating the “pulse on”$-$“pulse-off” analysis presented in Section 3.2.1 for the two relevant archival XMM-Newton observations. As before, here we can only make use of the EPIC-pn data from XMM-Newton for epochs X2 and X3. We fit the pulsed emission from all three epochs (X2, X3 and XN1) simultaneously, using the CUTOFFPL model again. Owing to the limited bandpass available for epochs X2 and X3, we link the high-energy cutoff between all epochs so that we can undertake a self-consistent comparison of the slope of the pulsed emission, given the degeneracy between $\Gamma$ and $E_{\rm{cut}}$ that can be seen in Figure 2. In Figure 7 we show the 2-D confidence contours for the photon index and the 2–10 keV flux for the pulsed emission from all three epochs. Although the pulsed flux does vary, in all three cases the spectral forms are clearly consistent. We therefore perform a final fit to the pulsed spectra in which $\Gamma$ is linked across all the datasets as well, leaving only the normalisation of the CUTOFFPL model free to vary between them. With this approach, we find $\Gamma=0.17^{+0.17}_{-0.19}$ and $E_{\rm{cut}}=4.7^{+0.9}_{-0.7}$ keV, consistent with the results obtained fitting the pulsed data from epoch XN1 only. The pulsed fluxes computed with this final analysis are given in Table 4.

4.2 Time-Averaged Analysis of the ULX State

We now undertake a simultaneous, multi-epoch analysis of the average spectra extracted from the epochs during which P13 exhibited ULX lumonisities. We take a similar approach to our analysis of the phase-resolved spectra from epoch XN1 (see Section 3.2.2) and fit the data from each epoch with three continuum components: a CUTOFFPL model to account for the pulsed emission, and an underlying 2-component continuum to account for the non-pulsed emission, again utilizing the DISKBB+DISKPBB combination. Having found that the spectral shape of the pulsed emission remains constant throughout all the epochs from which pulsations could be detected to date, we utilize this in our multi-epoch analysis and again fix the shape parameters ($\Gamma$, $E_{\rm{cut}}$) of the CUTOFFPL model to the best-fit values obtained from our analysis of the pulsed spectra, taking the values from our global fit to the pulsed emission from epochs X2, X3 and XN1 (see Section 4.1).

Given this approach, we focus here on epochs X2, X3 and XN1, since we are not able to resolve the pulsations in the Chandra data from epoch C1. However, we do note that the average spectrum from epoch C1 is very similar to epoch X3 (see Figure 5), so similar results would be expected. We fit the data from all three epochs simultaneously in order to investigate which of the spectral parameters drive the observed variability. This also allows us to make the simplifying assumption that the total flux from the accretion column is proportional to the pulsed flux. As such, we assume that the average CUTOFFPL flux in the 2–10 keV band (common to all epochs) shows the same relative variations as the pulsed flux from epoch-to-epoch (see Figure 7 and Table 4); the fluxes of the pulsed emission across all three epochs are therefore controlled by a single free parameter in our analysis, which we choose to be the pulsed flux in epoch XN1.

We also assume that the intrinsic neutral absorption column does not vary between the epochs considered (see e.g. Miller et al. 2013). Furthermore, during our analysis we find that the radial temperature index for the DISKPBB component ($p$) can be linked across all epochs; doing so only worsens the fit by $\Delta\chi^{2}=1$ for 2 fewer free parameters (an F-test444See Protassov et al. (2002) for appropriate uses of the F-test. confirms this improvement is not significant, with a false-alarm probability of 0.6). Both of the DISKBB and DISKPBB temperatures are found to vary. Linking either of these across all three epochs degrades the fit by $\Delta\chi^{2}>22$ (again for two fewer free parameters; here the F-test gives a false-alarm probability of $<3\times 10^{-5}$). The global fit for the model in which $N_{\rm{H,int}}$ and $p$ are linked across all epochs is very good, with $\chi^{2}$/DoF = 2149/2058, and the parameter constraints are presented in Table 5.

One notable result common to all our spectral analyses is that the radial temperature index for the DISKPBB component is rather steep; here we find $p>1.52$ (see also Sections 3.1 and 3.2.2). This implies that the emission from this component is significantly more peaked, i.e. dominated by a narrow range of temperatures, than a standard thin disk continuum (which should give $p=0.75$; Shakura & Sunyaev 1973). We therefore test here whether a range of temperatures is really required for this component, and replace DISKPBB with a single blackbody (BB) component. The fit in this case is actually slightly improved over the DISKPBB model, despite having one fewer free parameter: $\chi^{2}$/DoF = 2144/2059 (this is likely due to the upper limit of 2.0 we place on $p$ in our previous models). In contrast, the lower temperature component does require a range of temperatures; replacing the DISKBB component with a single BB worsens the fit by $\Delta\chi^{2}>30$ (for the same number of free parameters) regardless of the model used for the hotter component (DISKPBB or BB).

The results for the model in which the hotter component is treated as a single BB are also given in Table 5, and we show the data/model ratios for each epoch in Figure 8. We also show the relative contributions of the various components for the best-fit model in Figure 9; the evolution of the various components is qualitatively consistent with that expected from the evolution in the soft and hard band pulse fractions (Figure 6). Lastly, as a further sanity check, we compare the total flux for the CUTOFFPL model in epoch XN1 (against which the other epochs are scaled) obtained with this model to that expected based on the phase-resolved analysis of this epoch (Section 3.2.2) and the pulsed flux observed (Table 4). The flux for the CUTOFFPL component obtained for the minimum of the pulse cycle is $\sim$1.7 $\times 10^{-12}$ erg cm*-2* s*-1* in the 2–10 keV band (see Table 3). For a perfectly symmetric sinusoidal pulse profile, the total CUTOFFPL flux would be expected to be the sum of this flux and half of the pulsed flux, i.e. $\sim$2.6 $\times 10^{-12}$ erg cm-2 s*-1*. This is very similar to the total CUTOFFPL flux obtained for epoch XN1 in our multi-epoch spectral analysis with the DISKBB+BB+CUTOFFPL model (Table 5), particularly when accounting for statistical uncertainties and the fact that the pulse profile is not perfectly sinusoidal (Fürst et al. 2016b), suggesting that our spectral decomposition in this case is internally self-consistent.

4.3 The Low-Flux Observations

In addition to the ULX luminosity observations, we also undertake a simple analysis of the average spectra extracted from the observations taken during the off-state (epochs C2, C3, C4 and X1). Although there are some flux variations, the spectra from these epochs are all broadly consistent (key spectral parameters agree within their 90% errors when fit individually), so to maximise the S/N we co-add the Chandra data to form a single spectrum using ADDASCASPEC. The combined Chandra spectrum from epochs C2+C3+C4 and the EPIC-pn data from epoch X1 are shown in Figure 5. In order to get the best constraints for this state, we fit these Chandra and XMM-Newton datasets simultaneously with some simple models. We first apply a simple absorbed powerlaw continuum, fixing the intrinsic neutral absorption to $N_{\rm{H;int}}=8\times 10^{20}$ $\rm\thinspace cm^{-2}$ based on the fits to the broadband data, and assuming a common photon index for the Chandra and XMM-Newton datasets. However, we find that this does not fit the data ($\chi^{2}$/DoF = 63/29) and leaves obvious residuals. These are particularly prominent around $\sim$1 keV (see Figure 10), and are reminiscent of the contribution from thermal plasma emission at low energies.

We therefore add a MEKAL plasma component to the model. We assume this is only absorbed by Galactic absorption, that this component has solar abundances, and that the temperature remains constant between the Chandra and XMM-Newton datasets, but we also find that the normalisations for this component are consistent between the datasets, and therefore also link this parameter in our final model. The only parameter that is allowed to vary between the XMM-Newton and Chandra spectra is the normalisation of the powerlaw component. This model fits the off-state data very well ($\chi^{2}$/DoF = 28/27), and resolves the $\sim$1 keV residuals. We stress that both the MEKAL and the powerlaw components are required to fit the data; removing the latter results in a significantly worse fit again ($\chi^{2}$/DoF = 61/29). The best-fit model for the Chandra spectrum is shown in Figure 10, and the parameter constraints are given in Table 6.

The MEKAL component contributes at low-energies, fitting the $\sim$1 keV feature with iron L emission, but in addition to this a very hard powerlaw continuum ($\Gamma=0.7^{+0.4}_{-0.5}$) is required at higher energies, which we assume to be related to residual accretion onto the neutron star. However, the fact that the normalizations for the MEKAL component are consistent for the XMM-Newton and Chandra datasets, despite their significantly different extraction regions, suggests that this is not faint, extended emission from the broader galaxy, but is likely also be associated with P13. However, even if this is the case it is not clear that this emission is directly powered by accretion, so we compute the observed fluxes for both the total emission and just the powerlaw component (see Table 6). The latter provides the majority of the observed flux in the 0.3–10.0 keV band in both the Chandra and the XMM-Newton datasets; the residual accretion power is orders of magnitude lower than the peak flux observed from epoch XN1, regardless of the power source for the plasma emission.

5 The Broadband ULX Sample

For comparison with NGC 7793 P13, we also compiled the sample of ULXs with broadband spectra available in the literature at the time of writing. We focus only on sources for which NuSTAR has provided a robust detection out at least 20keV, which are not significantly confused with any other X-ray sources, and for which simultaneous data from at least one other mission with high S/N soft X-ray coverage (XMM-Newton and/or Suzaku) exist. This sample includes: Circinus ULX5 (Walton et al. 2013), Holmberg II X-1 (Walton et al. 2015b), Holmberg IX X-1 (Walton et al. 2014, 2017; Luangtip et al. 2016), IC 342 X-1 and X-2 (Rana et al. 2015), NGC 1313 X-1 (Bachetti et al. 2013; Miller et al. 2014; Walton et al. 2016b), NGC 5204 X-1 (Mukherjee et al. 2015) and NGC 5907 ULX (Walton et al. 2015a; Israel et al. 2017a; Fürst et al. 2017). Where necessary, we reprocessed the data with the latest calibrations, largely following the reduction procedures outlined in the referenced works, and also included the NuSTAR ‘spacecraft science’ data (see Section 2). Notable exclusions are NGC 1313 X-2, which is not detected much above 10 keV (Bachetti et al. 2013), and M82 X-1 and X-2, which are strongly blended in the NuSTAR data (e.g. Bachetti et al. 2014; Brightman et al. 2016a, b). We also do not consider NGC 5643 X-1 or M51 X-8; both are only weakly detected by NuSTAR above 10 keV in short observations, and neither of the NuSTAR observations of these sources have simultaneous high S/N soft X-ray coverage (Annuar et al. 2015; Krivonos & Sazonov 2016; Earnshaw et al. 2016). Finally, we also note that there is a further NuSTAR observation of IC 342 that we also do not include, as again this has no simultaneous high S/N soft X-ray coverage (Shidatsu et al. 2017). The majority of the sources considered here are also included in the ULX sample recently presented by Pintore et al. (2017).

The time-averaged broadband spectra for the selected sample, including NGC 7793 P13, are shown in Figure 11 (in the case of Holmberg IX X-1, we only show three of the six broadband observations available, corresponding to epochs 1, 2 and 5 as defined in Walton et al. 2017; the remaining three are similar to either the low or medium-flux states shown). There is clearly some diversity among the observed spectra; the lowest energies in particular are influenced by different levels of neutral absorption. However, there are also many similarities. For all the sources with low absorption, evidence for two thermal continuum components can be seen in the spectra below 10 keV (although the relative contributions of these components can vary both between sources and over time for an individual source, also helping to increase the diversity of the observed spectral shapes), similar to P13. Motch et al. (2014) also note that the spectrum of P13 below 10 keV is fairly representative of the broader ULX population. The qualitative similarity of the NuSTAR data in all cases is particularly striking; all the high-energy spectra break in a fairly narrow energy band above $\sim$3 keV to a steep continuum above 10 keV. Evidence for an additional high-energy powerlaw tail (or at least, an additional high-energy continuum component) has now been seen in several other systems in addition to P13 (Circinus ULX5: Walton et al. 2013; Holmberg II X-1: Walton et al. 2015b; Holmberg IX X-1: Walton et al. 2014, 2017; NGC 5204 X-1: Mukherjee et al. 2015; NGC 5907 ULX: Fürst et al. 2017). Despite the clear diversity among the spectra, the basic spectral components present in the broader ULX population, for which the nature of the accretors generally remains unknown, therefore appear similar to NGC 7793 P13 (see also Pintore et al. 2017).

6 Discussion

We have performed a detailed spectral analysis of the ultraluminous X-ray source NGC 7793 P13, focusing on the first high-energy detection of the source by NuSTAR. The recent detection of X-ray pulsations (Fürst et al. 2016b; Israel et al. 2017b) firmly identify it as a neutron star accretor, meaning the peak luminosity observed ($L_{\rm{X,peak}}\sim 10^{40}$ erg s*-1*) is highly super-Eddington. The low level of absorption towards P13 and the lack of source confusion makes it a key laboratory for understanding the super-Eddington accretion onto neutron stars, and also their relation to the broader ULX population. Through a combination of time-averaged, phase-resolved and multi-epoch analyses, we find that three spectral components are required to fit the data: two thermal blackbody components that contribute below $\sim$10 keV (with temperatures of the hotter and cooler components varying from $\sim$1–1.5 keV and from $\sim$0.3–0.5 keV between the observed epochs, respectively), and a third component that extends to higher energies. As discussed in Section 5, these components are qualitatively similar to those seen in other ULXs with broadband coverage.

6.1 The Thermal Components

Here, we examine the behaviour of these thermal components in order to try and determine their origin. One subtle difference between P13 and the broader ULX population worth noting is that for the hotter of the two components we find that the data prefer a steep radial temperature index (i.e. $p\gtrsim 1.0$) when fit with a disk model (DISKPBB). This is true for all our various analyses of the P13 data (the time-averaged and phase-resolved analysis of epoch XN1, and the multi-epoch analysis of epochs X2, X3 and XN1). In fact, for our multi-epoch analysis, the fit is slightly improved if we do not allow for a range of temperatures at all, and fit the hotter component as a single-temperature blackbody (Section 4.2). For other ULX systems, the radial temperature indices obtained for the hotter component have typically been much flatter, with $p<0.75$ (e.g. Walton et al. 2015b, 2017).

In Figure 12 we plot the flux of each of the DISKBB and BB components (computed over 0.01–100 keV, a sufficiently broad band to effectively be bolometric) against their temperatures. With only three data points to date, we do not formally fit any relations to these data, but for comparison we do overlay illustrative $F\propto T^{4}$ relations as expected for blackbody emission from a constant emitting area. The lower energy DISKBB component does appear to follow $F\propto T^{4}$ reasonably well, suggesting a roughly constant emitting area, but the same is not true for the BB component. In the latter case, epoch XN1 sits significantly off any extrapolation of an $F\propto T^{4}$ relation from the lower flux data, suggesting the emitting area for this component may be variable.

For the blackbody component, the maximum and minimum flux/temperature combinations (from epochs XN1 and X2, respectively) imply emitting areas with radii of $\sim$40–70 km from standard blackbody theory, or equivalently $\sim$20–35 $R_{\rm G}$ for a 1.4 $M_{\odot}$ neutron star (where $R_{\rm G}$ = $GM/c^{2}$ is the gravitational radius). This is comfortably larger than typical radii for the entire neutron star ($\sim$10 km), suggesting that this component must come from the accretion flow rather than from the surface of the neutron star itself, particularly given that the emission from the surface of a magnetised neutron star should primarily arise from small hotspots around the magnetic poles onto which the accretion is channelled. For the lower temperature DISKBB component, the fact that this component is required to have a range of temperatures already suggests that this must also be associated with the accretion flow. This is further supported by the emitting radius implied by the DISKBB normalizations, which are proportional to $R_{\rm{in}}^{2}\cos\theta/f_{\rm{col}}^{4}$. Here, $R_{\rm{in}}$ and $\theta$ are the inner radius and the inclination of the disc, and $f_{\rm{col}}$ is the colour correction factor relating the observed ‘colour’ temperature to the effective temperature at the midplane of the disc ($f_{\rm{col}}=T_{\rm{col}}/T_{\rm{eff}}$). This latter quantity provides a simple empirical correction accounting for the complex physics in the disk atmosphere. Even assuming that both $\cos\theta$ and $f_{\rm{col}}$ are unity (*i.e. *no colour correction and a face-on disc; note that this is also implicitly assumed for the hotter thermal component in our use of standard blackbody theory above), which would give the smallest radii, we find that $R_{\rm{in}}\sim 280$-360 km ($\sim$140–180 $R_{\rm G}$). This is at least an order of magnitude larger than typical neutron star radii. We will re-visit these radii in Section 6.3.

6.2 The Pulsed Emission/High-Energy Tail

In addition to the two thermal components, our analysis of the time-averaged spectrum finds that a third continuum component is required above 10 keV, as all two-component thermal models leave a clear excess at the highest energies probed by NuSTAR (Figure 1). As discussed previously, similar high-energy excesses have been seen in the average spectra of several other ULX systems now. In these previous works, the excess has typically been modelled as steep ($\Gamma\sim 3$) high-energy powerlaw tail to the lower-energy thermal continuum components, assumed to arise through Compton up-scattering similar to the X-ray coronae seen in sub-Eddington systems. When considering the time-averaged spectrum for P13, the same model fits the data well.

However, through phase-resolved analyses of the P13 data, we find that this excess is actually associated with the pulsed emission component arising from the magnetically channelled accretion columns. rather than from a traditional X-ray corona. The emission from the accretion column dominates over the two thermal components above 10 keV, naturally explaining the fact that the pulsed fraction increases with increasing energy below 10 keV, and appears to saturate at higher energies (Fürst et al. 2016b). This emission is well described with a simple CUTOFFPL model, which has a very hard rise ($\Gamma\sim 0.15$) below 10 keV before breaking to a steep spectrum above 10 keV thanks to a low cutoff energy ($E_{\rm{cut}}\sim 4.5$ keV). This is qualitatively similar to the pulsed emission from the the first ULX pulsar discovered, M82 X-2 (Brightman et al. 2016a, although for M82 X-2 the rise is not quite as hard and the cutoff at a higher energy, see Figure 2). The average spectrum of X-2 cannot easily be disentangled from M82 X-1, complicating spectral analyses. Nevertheless, this does show that the emission from the accretion column also extends up to these energies, and the pulsed fraction for M82 X-2 is also known to increase with increasing energy (Bachetti et al. 2014), despite the confusion with X-1. We therefore expect that the emission from the accretion column also dominates at the highest energies in this system.

6.3 Super-Eddington Accretion onto Pulsar ULXs

For super-Eddington accretion onto a compact object, the inner regions of the accretion flow are expected to have a large scale-height, owing to the increased support by the intense radiation pressure (e.g. Shakura & Sunyaev 1973; Abramowicz et al. 1988; Poutanen et al. 2007; Dotan & Shaviv 2011). This occurs interior to the point at which the accretion flow is locally at the Eddington limit, a radius commonly referred to as the ‘spherization radius’ ($R_{\rm{sp}}$). The geometry of a super-Eddington accretion flow is expected to resemble a funnel-like structure close to the central accretor that transitions to a more standard thin disk at larger radii. The central accretors powering ULXs have generally been assumed to be black holes in the literature, and theoretical effort has therefore primarily focused on super-Eddington accretion onto these objects, but King (2008) argue that the same framework should hold for super-Eddington accretion onto a neutron star (at least, in regions of the disk not significantly influenced by the neutron star’s magnetic field).

In the context of this model, the two thermal components required to successfully reproduce ULX spectra below 10 keV are often interpreted as arising from different regions in a super-Eddington accretion flow, with the hotter and cooler components arising from inside and outside the inner funnel, respectively (e.g. Middleton et al. 2015; Luangtip et al. 2016, but see Miller et al. 2014 for an alternative scenario). Large scale-height, optically-thick accretion flows are expected to show $p<0.75$ owing to the significant photon advection that should occur in this region (Abramowicz et al. 1988), as is often seen for the hotter component in ULX spectral modelling (e.g. Walton et al. 2015b, 2017), but is not seen for P13.

For ULX pulsars, the influence of the neutron star’s magnetic field also needs to be considered, This will disrupt the disk at the magnetospheric radius ($R_{\rm{M}}$), at which point the accreting material is forced to follow the field lines down onto the neutron star. Since the discovery of the first ULX pulsar (Bachetti et al. 2014), an extremely wide range of magnetic fields have been proposed in the literature in order to explain the observed properties of these sources, from $B\sim 10^{9-14}$ G (e.g. Kluźniak & Lasota 2015; Mushtukov et al. 2015). This naturally leads to a range of predictions for the accretion geometry, depending on the relative sizes of $R_{\rm{sp}}$ and $R_{\rm{M}}$. If $R_{\rm{M}}$ $>$ $R_{\rm{sp}}$ (high B-field), then the magnetic pressure will terminate the accretion flow before the large scale-height inner funnel can form. Alternatively, if $R_{\rm{M}}$ $\ll$ $R_{\rm{sp}}$ (low B-field), then the disk will only truncate very close to the neutron star, allowing the inner funnel to form unhindered, and the geometry will resemble that described above. Finally, for more intermediate configurations in which $R_{\rm{M}}$ $\lesssim$ $R_{\rm{sp}}$ (*i.e. *smaller, but still roughly similar), then the inner funnel should begin to form, but will be terminated prematurely by the magnetic field.

The fact that even after accounting for the pulsed emission from the accretion column we still require two thermal components to reproduce the P13 data supports the idea that there are two distinct regions to the accretion flow beyond $R_{\rm{M}}$, suggesting that $R_{\rm{M}}$ $<$ $R_{\rm{sp}}$. One might therefore expect to be able to associate the characteristic radii of the hotter and cooler components with $R_{\rm{M}}$ and $R_{\rm{sp}}$, respectively. Our previous crude estimates of these radii are qualitatively consistent with this scenario, with the hotter component arising from a smaller region. Interestingly, the radii of these two components appear to be roughly comparable (differing by a factor of $\sim$6; Section 6.1).

Several additional factors need to be considered to assess whether these radii could really be similar in reality. The values presented in Section 6.1 were estimated by setting both $f_{\rm{col}}$ and $\cos\theta$ to unity in order to assess whether either component could be associated with the neutron star surface. Both of these quantities can have a strong influence on the absolute radii inferred. However, as long as $f_{\rm{col}}$ and $\theta$ are relatively similar for both components, their relative radii should be fairly insensitive to the values adopted. On some level, we might expect that $f_{\rm{col}}$ should be larger for the hotter component (Shimura & Takahara 1995), which would increase the radius of this component relative to the cooler one. Large differences in $f_{\rm{col}}$ are not likely though (Shimura & Takahara 1995; Davis et al. 2006), so we expect this will only influence the relative sizes inferred by up to a factor of a few. In addition, if the BB component does arise from the regions of the flow interior to $R_{\rm{sp}}$, then this will likely be subject to some level of geometric collimation by the inner funnel. As such, the observed flux will be enhanced relative to the intrinsic flux by a beaming factor $b$ (i.e. $L_{\rm{int}}=bL_{\rm{obs}}$ adopting the nomenclature of King 2008, where $b\leq 1$, and $b=1$ represents isotropic emission). Since the outer disc would be largely immune to similar geometric beaming this would serve to decrease the radius of the hotter component relative to the cooler one, but only by a factor of $b^{-1/2}$. For illustration, the beaming factor of $b\sim 0.1$ suggested by King (2008) for neutron star accretors reaching $L_{\rm{X}}\sim 10^{40}$ erg s*-1* (similar to P13) would only change the relative radii of the thermal components by a factor of a few. Furthermore, truncation of the disc by the magnetic field will reduce the degree of geometric beaming relative to the case in which the magnetic field has a negligible effect on the flow (as assumed by King 2008). Given that these two issues effects act in opposition to each other, and may result in corrections of a similar order, the conclusion that $R_{\rm{sp}}$ and $R_{\rm{M}}$ are comparable would still appear to hold.

If $R_{\rm{M}}$ is fairly similar to $R_{\rm{sp}}$, then this could potentially explain why the hotter thermal component in P13 does not require a strong range of temperatures. If the magnetic field quickly truncates the disk for $R<$ $R_{\rm{sp}}$, then the geometrically thick inner region of the flow would only exist over a fairly narrow range of radii, and subsequently only emit over a narrow range of temperatures. The preference for a single blackbody/steep radial temperature index is likely driven in part by the assumption that the disc extends out to radii $\gg$ $R_{\rm in}$ (taken to be $R_{\rm{M}}$ in our simple picture), which is implicit in the DISKPBB model. This could artificially force $p$ to be steep in the fits since this is the only way to make the DISKPBB model effectively emit over a narrow range of temperatures, even if in reality $p<0.75$ between $R_{\rm{sp}}$ and $R_{\rm{M}}$, as would be expected for a thick super-Eddington disc.

It is interesting to note that King et al. (2017) also suggest that $R_{\rm{sp}}$ and $R_{\rm{M}}$ are similar for the three ULX pulsars identified to date. However, we caution that this is based on a model which assumes that the lower temperature component represents the blackbody emission from within the inner funnel of a super-Eddington flow. Furthermore, the model is designed to reproduce a luminosity–temperature relation of the form $L\propto T^{-4}$ (note the sign on the exponent), which has previously been claimed observationally for this component. This is primarily based on the work of Kajava & Poutanen (2009), who find an anti-correlation between luminosity and temperature with this approximate form by fitting the results for the soft thermal component compiled from a sample of ULXs. More recent analyses do not support such a relation for the soft component when individual sources are considered (Miller et al. 2013), instead finding positive indices (*i.e. *luminosity positively correlated with temperature). Critically, for P13 we find that both the cooler DISKBB and the hotter blackbody components show positive correlations between luminosity and temperature, and certainly do not follow $L\propto T^{-4}$. Nevertheless, our results would appear broadly consistent with the broad conclusions of King et al. (2017) with regard to the similarity of $R_{\rm{sp}}$ and $R_{\rm{M}}$.

Although it is challenging to constrain the absolute value of $R_{\rm{M}}$ for P13 from our analysis given the issues related to beaming and the colour correction discussed above, if $R_{\rm{M}}$ $<$ $R_{\rm{sp}}$ (as is supported by the requirement for two non-pulsing blackbody components) then we can obtain a rough upper limit on the magnetic field of the neutron star by taking the limit that $R_{\rm{M}}$ = $R_{\rm{sp}}$. For a dipole field and a typical neutron star mass and radius, the magnetospheric radius is given by $R_{\rm{M}}$ $\sim(2.6\times 10^{8})\dot{M}_{17}^{-2/7}B_{12}^{4/7}$ cm (Lamb et al. 1973; Cui 1997; Fuerst et al. 2017), where $\dot{M}_{17}$ is the accretion rate at $R_{\rm{M}}$ in units of $10^{17}$ g s*-1* and $B_{12}$ is the magnetic field strength in units of $10^{12}$ G. Standard super-Eddington theory suggests the location of $R_{\rm{sp}}$ is related to the input accretion rate: $R_{\rm{sp}}/R_{\rm{G}}\sim(27/4)\dot{m}_{0}$ (independent of the nature/magnetic configuration of the central accretor), where $\dot{m}_{0}$ is the input accretion rate in Eddington units (Shakura & Sunyaev 1973), so taking this limit allows us to estimate both $R_{\rm{M}}$ and $\dot{M}_{17}$, assuming conservatively that the accretion rate at both $R_{\rm{sp}}$ and $R_{\rm{M}}$ is the same as $\dot{m}_{0}$ (in reality, we expect there to be strong mass loss from the system; this should primarily occur interior to $R_{\rm{sp}}$, being the point at which the accretion flow becomes locally Eddington, and will reduce the accretion rate at $R_{\rm{M}}$ relative to that at $R_{\rm{sp}}$, which further emphasises that our calculation is an upper limit).

Our prior estimate of $R_{\rm{sp}}$ (Section 6.1) assumed both $\cos\theta$ and $f_{\rm{col}}=1$, to conservatively compare with neutron star radii. The former may be reasonable, since we appear to have an unobscured view of the inner regions of the flow (Sutton et al. 2013; Middleton et al. 2015), but the colour correction for a standard disc (expected beyond $R_{\rm{sp}}$) is more typically expected to be $f_{\rm{col}}\sim 1.7$ (Shimura & Takahara 1995). Adopting this value for $f_{\rm{col}}$ increases the estimate for $R_{\rm{sp}}$ to $\sim$900 km, or equivalently $R_{\rm{sp}}/R_{\rm{G}}\sim 450$, and therefore $\dot{m}_{0}\sim 66$. For a standard 1.4 $M_{\odot}$ neutron star the Eddington accretion rate is $\dot{M}_{\rm{E}}=2.2\times 10^{18}$ g s*-1* (assuming a radiative efficiency of $\eta\sim 0.09$), so the accretion rate at $R_{\rm{M}}$ would be $\dot{M}\sim 1.4\times 10^{20}$ g s*-1*. Substituting all this back into the original expression for $R_{\rm{M}}$ and rearranging, we therefore find that $B_{12}\lesssim 6$. Interestingly, the co-rotation radius for P13 is also $R_{\rm{co}}$ $\sim$ 900 km, so we have $R_{\rm{sp}}$ $\sim$ $R_{\rm{co}}$. Since we know that $R_{\rm{M}}$ $<$ $R_{\rm{co}}$, as accretion must be occurring in P13, this constraint would also imply $B_{12}\lesssim 6$. Although we expect the similarity of $R_{\rm{sp}}$ and $R_{\rm{co}}$ is coincidental (as long as $R_{\rm{M}}$ $<$ $R_{\rm{co}}$ there is no reason $R_{\rm{co}}$ should be an important radius in terms of the observed emission), this does provide further support for our conclusion that $R_{\rm{M}}$ $<$ $R_{\rm{sp}}$, and that the thick inner disc does form over some range of radii.

This is consistent with the estimate in Fürst et al. (2016b), who found $B_{12}\sim 1.5$ based on the observed spin-up rate of the neutron star. However, we stress that our calculation fundamentally depends on the association of the characteristic radius of the cooler thermal component with $R_{\rm{sp}}$. Furthermore, since $B\propto R_{\rm{M}}^{7/4}$, any errors in the radius will be magnified in the magnetic field limit estimated here. Our limit on $B$ therefore depends very strongly on the value assumed for the colour correction factor, since the radii of the thermal components are themselves proportional to $f_{\rm{col}}^{2}$. We also stress that $R_{\rm{M}}$ is most sensitive to the dipolar field, as any quadrupolar component will decay much more rapidly with radius, so these constraints should be considered on the dipolar component only. A stronger quadrupolar component to the magnetic field that only acts close to the neutron star, as suggested by Israel et al. (2017a), cannot be excluded by these considerations.

One of the basic predictions of super-Eddington accretion is that powerful winds should be launched from the radiation-pressure dominated regions of the flow (*i.e. *primarily interior to $R_{\rm{sp}}$; e.g. Poutanen et al. 2007). Indeed, unambiguous evidence for powerful outflows has now been seen through resolved, blueshifted atomic features in a handful of other ULX systems with deep XMM-Newton coverage, primarily in the soft X-ray band ($<$2 keV) but in the case of NGC 1313 X-1 also in the iron $\rm{K}\alpha$ band (Pinto et al. 2016, 2017; Walton et al. 2016b). If our proposed scenario for P13 is correct, and thick inner regions of the flow do begin to form, then we might expect to see evidence for a similar wind in P13, either through ionised absorption features if the outflow crosses our line of sight, or ionised emission otherwise. Given the low levels of neutral absorption here, this could be tested with future deep observations combining either XMM-Newton or XARM (the Hitomi replacement mission; Takahashi et al. 2012) together with NuSTAR.

6.4 The Low-Flux State

In addition to the bright state observations, we also investigate the spectrum obtained during the off-state seen between 2011 and 2013 (epochs C2, C3, C4 and X1). The spectrum during these epochs is still very hard, and remains variable (Figure 5). However, simple powerlaw fits leave clear line-like residuals at $\sim$1 keV, and the fit is significantly improved with the addition of a plasma component, which we argue is associated with the P13 system, rather than galactic-scale plasma emission (e.g. star formation). If P13 does launch a strong wind while at ULX luminosities, this plasma component could potentially arise from collisionally-ionised material previously ejected from the accretion flow that becomes visible with the extreme drop in flux from the central accretor.

The hard ($\Gamma\sim 0.7$) and variable powerlaw emission almost certainly suggests that there is still some residual accretion onto the neutron star during these off-states. After the discovery of pulsations from M82 X-2, Tsygankov et al. (2016) suggested that the low-flux states seen in this source could be related to the onset of the propeller effect. This occurs when $R_{\rm{M}}$ moves outside $R_{\rm{co}}$, such that the magnetic field ends up acting as a barrier to further accretion, rather than channelling the material onto the neutron star poles. As such, for $R_{\rm{M}}$ $>$ $R_{\rm{co}}$ the accretion rate, and therefore the observed X-ray luminosity, suddenly drops precipitously. We do not see any clear evidence for an accretion disk component in these data; if we replace the MEKAL component with DISKBB, the fit is better than the simple powerlaw model ($\chi^{2}$/DoF = 44/27), but not as good as the MEKAL model, and the line-like residuals at $\sim$1 keV remain. This could suggest the disk would need to be truncated at a sufficiently large radius that it makes a negligible contribution to the observed flux in the X-ray band.

Nevertheless, the residual accretion flux is still quite substantial, with $L_{\rm{X}}\sim 3-6\times 10^{37}$ erg s*-1* in the 0.3–10.0 keV band, i.e. $\sim$0.15–0.3 $L_{\rm{E}}$ for a 1.4 $M_{\odot}$ neutron star. This is quite similar to M 82 X-2, which still radiates at $\sim$1038 erg s*-1* in its off-state (Brightman et al. 2016a). As noted by Tsygankov et al. (2016), this only requires a few percent of the accreting matter to leak through the magnetosphere, which may not be unreasonable. However, there is now evidence that the low-flux periods in M82 X-2 are related to its $\sim$62 d super-orbital cycle (Brightman et al. 2017, submitted), which would call into question the propeller interpretation for that system. If confirmed, this would establish a significant difference between M82 X-2 and both NGC 5907 ULX and P13, as the timing and duration of the off-states seen in the latter two means they cannot be related to their similar $\sim$64 d and $\sim$78 d periods. This means that there may not be a single physical origin for these off-states common to all three systems, despite the changes in observed flux being similar in all three cases, and so the propeller effect cannot be excluded for the low-flux states in P13.

6.5 Implications for the ULX Population

The qualitative similarity of the broadband spectrum of P13 to the rest of the ULX population has significant implications. Since P13 is a known super-Eddington accretor, the XMM-Newton+NuSTAR observation of this source presented here robustly demonstrates that these spectra really are associated with super-Eddington accretion. This confirms the conclusions of our previous work on broadband observations of ULXs (Bachetti et al. 2013; Walton et al. 2014, 2015b, 2017; Rana et al. 2015; Mukherjee et al. 2015), which were based on the lack of similarity between ULX spectra and the standard modes of sub-Eddington accretion onto black holes (see Remillard & McClintock 2006 for a review), rather than a well-established knowledge of the spectral appearance of super-Eddington accretion at ULX luminosities.

The similarity between P13 and the broader ULX population may even suggest that more of these sources could be powered by accreting neutron stars. Although pulsations have not currently been observed from any other sources (Doroshenko et al. 2015), for M82 X-2 the pulsations are known to be transient (Bachetti et al. 2014), and this also appears to be the case for NGC 5907 ULX (Israel et al. 2017a). Pintore et al. (2017) also suggest that additional ULXs could be pulsar systems, arguing that ULX spectra are well fit by spectral models typically applied to sub-Eddington pulsars: a low-energy blackbody (peaking below 2 keV) and a powerlaw with an exponential cutoff produced by an accretion column that dominates at higher energies ($>$2 keV). Our work suggests that the spectral decomposition used by Pintore et al. (2017) is not actually correct, and that the situation in P13 is more complex as the 2–10 keV emission does not exclusively arise from the column. Instead, we find that this only reliably dominates the emission above 10 keV, resulting in the high-energy excess seen with models that only invoke two continuum components (Figure 1). It is possible therefore that the presence of such a hard excess is indicative that the accretion in the innermost regions occurs via magnetically channelled columns, rather than via a disc. We stress again that a number of other ULXs have now been seen to exhibit similar excesses above $\sim$10 keV, so if this is correct then neutron stars/pulsars could potentially be prevalent in the ULX population.

Since the pulsations are seen to be transient in 2/3 of the known ULX pulsar systems, additional means of firmly identifying pulsar/neutron star accretors may well be the key to expanding this population. Indeed, King et al. (2017) suggest that pulsations may only be observable from neutron star ULXs when $R_{\rm{sp}}$ and $R_{\rm{M}}$ are similar (as we find appears to be the case for P13), as they argue that for $R_{\rm{M}}$ $\ll$ $R_{\rm{sp}}$ the pulsations could be diluted to the point of being undetectable by the non-pulsed emission from the disc. Although our work finds that the highest energies ($\gtrsim$15 keV) should be relatively immune to such issues as the column always dominates, the low count-rates observed at these energies would still serve as a major barrier to the detection of pulsations, further emphasising the need for additional means of identification.

One alternative method may be to search for other sources that show similar off-states to those seen from all three ULX pulsars, as suggested by Fürst et al. (2016a). In addition to this, based on our analyis of P13 in this work we tentatively suggest that finding that the hotter thermal component does not require a run of temperatures may be another indicator of a neutron star accretor. This would provide an indication that the disc is truncated close to $R_{\rm{sp}}$, despite the source maintaining ULX luminosities, and would in turn imply a magnetized central accretor.

7 Summary and Conclusions

We have undertaken a detailed X-ray spectral analysis of the available data for P13 – a key laboratory for ULX pulsars owing to its lack of source confusion, persistently detected pulsations and low obscuration – focusing on the first high-energy detection of this source with NuSTAR. Through time-averaged, phase-resolved and multi-epoch studies, we find that two thermal components are required to fit the data below 10 keV, in addition to a third continuum component above 10 keV to account for the high-energy excess seen with dual-thermal models. This is qualitatively very similar to the rest of the ULX population with broadband coverage, for which the nature of the accretors generally remain unknown, confirming the super-Eddington nature of these sources and suggesting that additional ULXs may be neutron stars/pulsars. Our phase-resolved analysis finds that the third, high-energy component is associated with the magnetically collimated accretion columns that must be present in P13 in order to produce the observed pulsations. Neither of the two thermal components can arise from the neutron star itself, as their characteristic radii at least an order of magnitude too large, suggesting that they both arise from the accretion flow. In turn, this suggests that the accretion flow has two distinct regions, as expected for super-Eddington accretion flows which should transition from a standard thin outer disc to a geometrically thick inner disc at the point that the flow reaches the local Eddington limit ($R_{\rm{sp}}$). However, the radii for these two components appear to be comparable, and the hotter component is better described as a single blackbody than a thick, advection-dominated accretion disc, which should have a strong run of temperatures.

We suggest that P13 may represent a scenario in which the magnetospheric radius – the point at which the disk is truncated and the accreting material is forced along the magnetic field lines – is smaller than, but still similar to $R_{\rm{sp}}$. This would allow the thicker inner regions of the super-Eddington flow to begin to develop, before being truncated prematurely by the magnetic field. Such a scenario would simultaneously explain the need for two thermal components, their fairly similar radii, and the fact that the hotter component does not need a run of temperatures, as the thick inner disk would only persist over a fairly narrow range of radii, and thus only emit over a fairly narrow range of temperatures. This similarity would imply an upper limit to the neutron stars dipolar magnetic field of $B\lesssim 6\times 10^{12}$ G, albeit with large uncertainties. If this is correct, evidence of similar truncation of the disc in other sources could offer a means of identifying additional neutron star/pulsar ULXs, as this would not be expected for black hole systems.

Finally, P13 is known to exhibit unusual ‘off’ states in which the X-ray flux drops by many orders of magnitude, potentially related to the propeller effect. We also examine the spectra obtained during one of these events, and find that the data require both a hard powerlaw component, suggesting residual accretion onto the neutron star, and emission from a thermal plasma, which we argue is likely associated with the P13 system.

ACKNOWLEDGEMENTS

The authors would like to thank the reviewer for the feedback provided, which helped to improve the final version of the manuscript. DJW and MJM acknowledge support from an STFC Ernest Rutherford fellowship, ACF acknowledges support from ERC Advanced Grant 340442, and DB acknowledges financial support from the French Space Agency (CNES). This research has made use of data obtained with NuSTAR, a project led by Caltech, funded by NASA and managed by NASA/JPL, and has utilized the NUSTARDAS software package, jointly developed by the ASDC (Italy) and Caltech (USA). This work has also made use of data obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States.