© ESO 2018

1. Introduction

Ultra-luminous X-ray pulsars (ULPs) are accreting neutron stars that show clearly detected pulsations and have inferred isotropic luminosities above L > 10³⁹ erg s⁻¹. These very high luminosities make them clearly super-Eddington, and therefore a challenge to our understanding of accretion. Currently, four ULPs are known, all located in nearby galaxies: M82 X-2 (Bachetti et al. 2014), NGC5907 ULX (Israel et al. 2017a), NGC300 ULX1 (also known as SN2010da; Carpano et al. 2018), and NGC7793 P13 (hereafter P13; Fürst et al. 2016; Israel et al. 2017b). Brightman et al. (2018) also identified the ultra-luminous source M51 ULX 8 as a likely neutron star accretor through the discovery of a cyclotron resonant scattering feature. However, no pulsations have been discovered from this source to date.

Before the discovery of pulsations, it was generally assumed that the compact object in ultraluminous X-ray sources (ULXs) was a black hole because the observed luminosities are more readily explained with a more massive compact object. To reach these luminosities in a sub-Eddington disk accretion regime, masses of 10^2–5 M_⊙ were implied (e.g., Miller et al. 2004). However, spectroscopic evidence, first by XMM-Newton (e.g., Stobbart et al. 2006; Gladstone et al. 2009) and later in particular by NuSTAR, has shown that for most ULXs stellar-remnant compact objects accreting at super-Eddington luminosities are more likely (e.g., Bachetti et al. 2013; Walton et al. 2013, 2014, 2015a,b; Rana et al. 2015; Mukherjee et al. 2015; Fürst et al. 2017; Walton et al. 2018a,b). For NGC7793 P13, Motch et al. (2014) could place a stringent constraint on the dynamical mass of the accretor of <15 M_⊙, confirming the super-Eddington explanation (of course, the later discovery of pulsations from P13 limits the maximum mass even further). The spectral similarity between ULPs and bright ULXs opens the possibility that most known ULXs are actually powered by an accreting neutron star, and that the expected pulsations are suppressed or diluted through a yet to be identified process (Pintore et al. 2017; Walton et al. 2018a).

Due to their large distances and sometimes high extinction, the optical counterparts of ULXs are difficult to identify. P13 is an exception, initially discovered in the X-rays by Read & Pietsch (1999), and whose mass donor has been identified as a B9Ia super-giant (Motch et al. 2011). Motch et al. (2014) identified an optical and UV photometric period of ≈64 d, which is also present in the radial velocity of the He II emission. While the origin of the He II emission line is debated (Fabrika et al. 2015), Motch et al. (2014) interpret the clearly detected period as the orbital period of the system and find a dynamical mass constraint of < 15 M_⊙.

A similar but significantly different period was identified in X-rays by Hu et al. (2017) at 65.05 ± 0.1 d, using long-term Swift/XRT monitoring data. Hu et al. (2017) discuss the difference between the optical/UV and X-ray period, and propose that it might be due to a beat frequency with a super-orbital period. The super-orbital period could be caused by a warped and precessing accretion disk (Hu et al. 2017). The arrival times of the optical maxima have been observed to vary, indicating the presence of a very long (∼2500 d) super-orbital period (Motch et al. 2014; Hu et al. 2017).

On the other hand, Hu et al. (2017) also discuss the possibility that both ∼65 d periods could actually be super-orbital in nature, and that the true orbital period is much shorter. This argument is based on the fact that the orbital periods for two other known ULPs are only on the order of a few days: 2.5 d for M82 X-2 (Bachetti et al. 2014) and likely ∼5 d for NGC5907 ULX (although in this case periods of up to a few months cannot be ruled out at the 3-σ level; Israel et al. 2017a). The observed phase jitter in the optical period in P13 could furthermore indicate an only semi-periodic behavior, which is expected for super-orbital periods but unlikely for orbital periods (Fürst et al. 2016). Middleton et al. (2018) suggest that this period is due to Lense–Thirring precession of the accretion flow, which might provide constraints on the equation of state of the neutron star.

Identifying the correct orbital period is an important step in understanding these systems, in particular with respect to formation and evolution models (Marchant et al. 2017). The best way to directly determine the orbital ephemeris in the X-ray band is to use the variation of the pulse period as a function of orbital phase. For geometries where we do not view the system face-on, the pulse period will vary periodically due to the Doppler effect.

With a pulse period of P ≈ 415 ms, P13 shows the fastest spin of all known ULPs, making it an ideal candidate for this search. The spin-period can be determined to within 10 ns with XMM-Newton with exposure times around 50 ks. The secular spin up Ṗ ≈ 3.5 ≈ 10⁻¹¹ s s⁻¹ has been well established over a three-year baseline (2013–2016, Fürst et al. 2016; Israel et al. 2017a).

To measure the influence of the orbital motion on P and Ṗ, and thereby to determine the ephemeris, we obtained five XMM-Newton and three NuSTAR observations, spread out semi-regularly over 65 d. This allowed us to determine the pulse period at different phases. These observations were supported by further XMM-Newton and NuSTAR observations, which increased the baseline to determine the secular spin-up trend. Table 1 gives an overview of the data used.

Here we present the timing analysis of these data, while the spectral analysis will be presented in a forthcoming work. In Sect. 2 we describe the observations and our data reduction. In Sect. 3 we briefly describe the analysis method, present the timing results, and describe them with a model for the orbital variation. In Sect. 4 we discuss the implications of our results. Uncertainties are given at the 90% confidence level unless otherwise noted.

2. Observations and data reduction

2.1. Swift

NGC7793 P13 has been frequently monitored by the Neil Gehrels Swift Observatory (Swift; Gehrels et al. 2004) across several periods, first for about ∼10 weeks in 2010, then for ∼6 months spanning late 2012 and early 2013, for another ∼5 weeks spanning late 2014 and early 2015, and then consistently since May 2016. Although there are occasionally longer gaps even within these monitored windows, the typical observing cadence is 5–10 days, and the average exposure is ∼2 ks. Figure 1 shows the long-term light curve of these observations obtained with the XRT and UVOT instruments (Burrows et al. 2005; Roming et al. 2005).

We extracted all available XRT data between 16 Aug 2010 (ObsID 00031791001) and 25 Jan 2018 (ObsID 00093149030) with the standard Swift/XRT processing pipeline (Evans et al. 2009), adding over 50 more observations to the light curve compared to the data presented by Hu et al. (2017). The data are binned such that there is a single 0.3–10 keV flux measurement from each observation. Aside from the 2012–2013 data, during which the source was in a low-flux state and was not detected within any individual observation (Motch et al. 2014; Fürst et al. 2016), an indication of a periodicity on the order of ∼65 days is visibly present in the X-ray data, with the observed count rate varying by a factor of ∼3–4 from peak to trough.

To extract the UVOT (Roming et al. 2005) data we downloaded all available Swift data of NGC 7793. We used a circular source region with a radius of 5″, as recommend by the Swift UVOT Software Guide¹ centered on the \chandra coordinates for P13 (RA = 23 h 57′ 50.9″, Dec = −32° 37′ 26.6″; Pannuti et al. 2011). As background region we chose a circular region with a radius of 15″, located 35″ east of P13 in a source-free region. As the source is located in the outskirts of the galaxy, some contamination of the magnitudes from the galactic background light is present. However, as we are only studying the variability of the source here, this should not influence our results. The UVOT data were extracted using the standard software as distributed with HEASOFT v6.20, in particular we used uvotimsum to sum up all exposures in each ObsID and uvotsource to extract the source magnitude.

2.2. NuSTAR

Data from NuSTAR (Harrison et al. 2013) were reduced using the standard pipeline, nupipeline, provided in the NuSTAR Data Analysis Software (v1.8.0), with standard filtering and NuSTAR CALDB v20171204. Source products were extracted from circular regions of radius 70″ for both focal plane modules (FPMA/B) using nuprodcuts, with the background measured from circular regions with a radius of 120″ on the same detectors as P13. Light curves were extracted in the 3–78 keV energy band with a maximum time resolution of 0.1 s. All time information was transferred to the solar barycenter using the DE-200 solar system ephemeris (Standish et al. 1992).

2.3. XMM-Newton

Data from XMM-Newton (Jansen et al. 2001) were reduced with the XMM-Newton Science Analysis System (v15.0.0), following the standard prescription.² Owing to its superior time resolution of 73.4 ms, in this work we only consider data from the EPIC-pn detector (Strüder et al. 2001). The data were taken in full frame mode and raw data files were cleaned and calibrated using epchain. Source products for the 2016 data were extracted from circular regions of radius ≈40″, following the same method as described in Fürst et al. (2016). For epochs 2017A–2017G we used a smaller circular region with a radius of 20″ to avoid contamination from another transient source close to P13. For the latest epochs we used a circular region with a radius of 30″ as the other source had faded. Background flaring was only problematic in epoch 2017A; for all other observations the full exposure time could be used (see Sect. 3.3). See Table 1 for a complete observation log.

3. Analysis

3.1. Optical and UV

We first performed a timing analysis of the UVOT data to confirm and update the optical period discovered by Motch et al. (2014), who mainly used data from ground-based telescopes. These authors found an optical period between ∼63.5 and 65.3 d, depending on the exact dataset used, with considerable phase jitter for the time of maximum light. Their best-fit solution implies a period of 63.52 d, with a super-orbital period of 2620 d, on which the phases of maximum light move around the average. Hu et al. (2017) used Swift/UVOT data up to December 2016 and found a period of 64.24 ± 0.13 d.

Most UVOT observations were performed with the “u” filter (central wavelength 3465 Å, 114 observations), but some were also performed using the “uvw1” (2600 Å, 19 observations), “uvw2” (1928 Å, 20 observations), and “uvm2” (2246 Å, 27 observations) filters; see Poole et al. (2008) for details about the UVOT photometric filters. The observed magnitudes of all the filters are comparable; only the uvm2 filter measures slightly lower values, as shown in the top panel of Fig. 1. Nonetheless, we decided to use all data together to search for periodicities to obtain the best signal statistically. Using only the u filter data gives consistent results, albeit with larger uncertainties.

We search for periodicities using a Lomb–Scargle periodogram (Scargle 1982) and epoch-folding (Leahy et al. 1983). To determine the uncertainties on the period, we simulate 5000 light curves based on the folded profile at the strongest period, with an added white noise term to obtain the same variance as in the original data (Davies 1990). For each simulated light curve we perform the same period search as for the real data, measure the strongest period, and report the 1-σ standard deviation of the distribution of these periods.

We find a strong signal at 63.9 ± 0.5 d, in very close agreement with the results by Motch et al. (2014) and Hu et al. (2017), as shown in the top panel Fig. 2. We fold the UVOT data on that period, and mark the times of maximum light based on the folded profile in Fig. 1. As can be seen, in the data after April 2016 (t = 2100 d) the times of maximum light are difficult to identify in the UVOT light curve. We therefore refrain from trying to fit a jitter of arrival times with a super-orbital period, as done by Motch et al. (2014) and Hu et al. (2017).

3.2. X-rays

As shown in the bottom panel of Fig. 1, the X-rays also show strong variability on a timescale similar to that of the optical data. Hu et al. (2017) identified an X-ray period of 65.05 ± 0.1 d after subtracting the average count rate of each observational epoch. This method of subtracting the average count rate is used to make the fainter observations in 2010 more comparable to the observations in 2015–2017. To measure the periodicity, we follow a slightly different approach here, and remove a linear brightening trend from the data. This approach is motivated by the fact that the average flux of P13 seems to be continuously brightening in the X-rays since the beginning of the observations. We also tried using a quadratic trend, but found no significant differences in the timing result.

Similar to our search in the UVOT data, we calculated the Lomb–Scargle periodogram and the epoch folding power using the detrended XRT light curve. The results are shown in the bottom panel of Fig. 2 and show a strong peak at 66.8±0.4 d (1-σ uncertainty). This period is significantly longer than the value given by Hu et al. (2017). This difference is due to our different approach of subtracting a linear trend from the data, and to the larger dataset used in this work. We then fold the XRT light curve on this period and mark the times of flux maximum in Fig. 1 based on the resulting profile. We note that this way we identify the data taken in December 2014 (t = 1600) as the falling flank of the X-ray variability, i.e., the time of maximum flux occurred just before the XRT observation took place, while Hu et al. (2017) assumed that the maximum flux occurred during the XRT observations.

The peak in the Lomb–Scargle periodogram and the epoch folding result is relatively broad, indicating that the X-ray period might not be perfectly stable. This might indicate that the X-ray period is actually only quasi-periodic or that it varies on a super-orbital timescale (see, e.g., Middleton et al. 2018). However, for the sake of simplicity we refer to it as the X-ray period in this paper.

Fig. 1. Panel a: light curve from Swift/UVOT, with one data point per observation. Black circles are u-filter, red squares are uvw1-filter, green triangles are uvw2-filter, and blue diamonds are uvm2-filter. Panel b: Swift/XRT light curve in the 0.3–10 keV energy band, with the same binning as the UVOT light curve. Dotted lines mark the times of maximum light in the UVOT data, assuming a constant period of Popt = 63:9 d and dashed lines mark the maximum of the folded XRT profile with a constant period of PX = 66:9 d.

Open with DEXTER

Fig. 2. Panel a: epoch-folding (black, left y-axis) and Lomb–Scargle Periodogram (green, right y-axis) of all Swift/UVOT data. Panel b: same as panel a but for the Swift/XRT data. The dotted vertical line marks the strongest period of the UVOT data (which is consistent with the implied orbital period from the pulse period analysis), the dashed line the period of the XRT data. The dashed horizontal line marks the 99.9% false alarm probability of the epoch folding results based on an F-distribution corrected for number of trial periods (Davies 1990).

Open with DEXTER

3.3. Pulsations

To determine the pulse period P we used all available XMM-Newton and NuSTAR data (see Table 1), apart from the XMM-Newton observations in 2012, which were taken during a very low flux state (Fürst et al. 2016). All Swift/XRT data have exposure times that are too short to identify the pulse period given the low average count rate of P13. For XMM-Newton we used EPIC-pn barycentered light curves extracted on the fastest time resolution of the detector of 73.4 ms and between energies of 0.3–10 keV. For NuSTAR we used barycentered and source-filtered event files (time-tagged with an intrinsic time resolution of 2 μs), with energies between 3 and 79 keV. To increase the signal-to-noise ratio (S/N), we combined the events of FPMA and FPMB.

As a first step, we calculated a simple power spectrum based on the fast Fourier transform for each dataset. However, no coherent pulsations were significantly detected in any observation in 2017 using these power spectra. This is contrary to the discovery of the pulsations, which showed up very clearly in the XMM-Newton data of 2012, 2013, and 2014 (Fürst et al. 2016; Israel et al. 2017b). The non-detections are probably due to the overall lower pulsed fractions and shorter exposure times of the new observations compared to previous observations.

In the next step we applied the epoch-folding technique, which is more sensitive to non-sinusoidal profiles. To narrow down our search range, we estimated the pulse period by extrapolating the secular spin-up as measured in our previous XMM-Newton and NuSTAR observations (Ṗ_sec ≈ 3.5 × 10⁻¹¹ s s⁻¹; Fürst et al. 2016). We searched around that period in a range of ±0.1 ms. In all observations besides 2017A this led to a clear detection of the pulse period.

To refine the determination of P and its derivative Ṗ, we then use an accelerated epoch-folding search, i.e., searching a 2D-grid in P-Ṗ space. We typically used a grid size of 200 × 200 points, centered on the best period found in the standard epoch folding, and on zero for Ṗ, with a search range of ±0.05 ms and ±10⁻⁹ s s⁻¹, respectively. As in Fürst et al. (2016), we determine the uncertainties on P and Ṗ from the full-width half-maximum (FWHM) contour of the χ² distribution. We show a typical example of this search in Fig. 3. We detect a period in each observation with a confidence larger than 99.99% (apart from epoch 2017A, see below). We calculated the significance based on a χ² distribution with 11 degrees of freedom (i.e., the number of phase bins in each pulse profile minus one, Davies 1990) and took the number of trial periods (200 in P and Ṗ each, i.e., a total of 40 000) into account.

Table 1 gives the measured P and Ṗ values in all observations. The first four measurements were published by Fürst et al. (2016) and the values for epochs X2/2013 and X3/2014 are also in good agreement with those found by Israel et al. (2017a). For the campaigns during which XMM-Newton and NuSTAR observed quasi-simultaneously (epochs 2016, 2017B, and 2017I), we performed independent searches in XMM-Newton and NuSTAR. While the values are fully consistent with each other, the NuSTAR observations provide the more precise measurement due to their longer duration. In the rest of the analysis we only use the NuSTAR value for those epochs.

Observation 2017A did not lead to a significant detection (at the 1-σ limit) of the pulse period. However, we find a peak in the epoch-folding result around 415 ms when selecting a window from between 5 and 15 ks after the start of the observation. The rest of the observation is plagued by strong background flaring, which likely dilutes the pulsed signal. While this detection is still not significant by itself, we report it as well in Table 1, but do not use it in our further analysis.

	Fig. 3. Example of the epoch-folding in P–Ṗ-space for NuSTAR observation 2017F. The χ2 value at each P–Ṗ pair is color-coded. The graphs at the top and on the right show the cut along the respective dimension for the most significant measurement at P = 415.7050 ms and Ṗ = 6.5 × 10−11s s−1. The orange dashed line shows the 99.99% significance line taking the number of trial periods into account.
Open with DEXTER

3.4. Pulse profiles

To investigate the behavior of the pulse profile and the pulsed fraction as function of time, we folded all available data on their respective pulse period (Fig. 4). The data were binned into 12 phase bins. We calculated the pulsed fraction PF as where p_i is the pulse profile of the respective observation.

As can be clearly seen in Fig. 4, the pulsed fraction during our 2017 campaign is on average much smaller than during the first three observations presented by Fürst et al. (2016). Furthermore, the pulse profiles seem to be less sinusoidal, even when the pulsed fraction is relatively high, i.e., in epoch 2017B which shows a narrower peak in both XMM-Newton and NuSTAR.

Fig. 4. Pulse profiles of all available XMM-Newton and NuSTAR observations, folded on the respective measured pulse period. The profiles are repeated once for clarity and phase-aligned by hand to have the minimum at phase 0. For epochs where XMM-Newton and NuSTAR observations where taken simultaneously, NuSTAR is shown in blue using the right-hand y-scale. Each panel is labeled with the epoch name, pulse period, and measured pulsed fraction in percent. For details see text.

Open with DEXTER

3.5. Determination of the orbit

As can be seen in Table 1 and Fig. 5b, the pulse periods around MJD 57900 show an apparent spin-down (between epochs 2017D to 2017G). This spin-down is a clear indicator of a strong Doppler effect due to the orbital motion, as a torque reversal on these timescales is unlikely. We therefore fit an orbital model to the data between 2016 and 2017. We note that we do not use data from observation 2017A, as the period could not be independently found there and only use the NuSTAR values from epochs 2016, 2017B, and 2017I given their much higher precision compared to XMM-Newton. We initially do not use older data, as over that long stretch of time a variable accretion torque is likely.

The free parameters in the model are the orbital period P_orb, the intrinsic pulse period P_spin, at a pre-defined time T₀ (which we fix at MJD 57530.00), the projected semi-major axis a sin i, the eccentricity e, the time of periastron τ, and the argument of periastron ω. Additionally, we allow for an ad hoc, constant spin-up parameter ξ. This model therefore has seven free parameters, compared to nine pulse period measurements, i.e., we only have two degrees of freedom. We do not find a good description of the data in terms of χ² (χ² = 4.9 for 2 d.o.f.); however, the secular spin-up trend and the indication of an orbital Doppler shift in 2017 is well captured. We find an orbital period around 63.7 d and a very small eccentricity (e ≈ 0.05). All parameters are given in Table 2.

The observed pulse periods will also be influenced by the transfer of angular momentum onto the neutron star from the accreted matter. Following the theory by Ghosh & Lamb (1979a,b), the observed pulsed period change is proportional to P_spinL^6/7 in the disk accretion case, where L is the intrinsic luminosity of the system. The exact proportionality factor depends on different factors, like the magnetic field strength and the way the accretion disk couples to the magnetosphere. For simplicity, we describe it with a simple spin-up parameter b, for which we can fit directly (Marcu-Cheatham et al. 2015; Bissinger 2016). The theory is based on a subcritical, geometrically thin accretion disk, and we therefore caution that the coupling might be significantly different in the super-Eddington case. Nonetheless, this approach should give us a good indication of the transfer of angular momentum in the system.

We use the data from the Swift/XRT monitoring as a tracer of the intrinsic luminosity (Fig. 5a). As the Swift/XRT monitoring has relatively large gaps, however, we do not use the data directly, but interpolate the light curve folded on the 66.9 d X-ray period plus a linear brightening trend (see Sect. 3.2). This interpolation removes short-term variability, but provides a good description of the overall behavior of the X-ray light curve (see Fig. 5). We thereby infer that the X-rays behaved in a similar way during the gaps of the XRT monitoring. As the spin-up depends on the luminosity, we give the proportionality factor b in units of Ṗ at an XRT count rate of 0.08 cts s⁻¹, which corresponds roughly to a luminosity of L₀ = 3.8 × 10³⁹ erg cm⁻² s⁻¹ (for a distance of 3.6 Mpc).

With this approach we find a very good description of the pulse period evolution between 2016–2017 (Fig. 5). The best-fit orbital period is determined to be d, and again the eccentricity is relatively small, . All parameters are summarized in Table 2.

We find that the orbital period is very close to the period we found in the optical and UV data. This indicates that the optical photometry indeed varies mainly due to orbital effects, as discussed by Motch et al. (2014). These authors assume that the variability is due to the X-ray illuminated side of the optical companion rotating in and out of view.

On the other hand, the timing solution clearly shows that the X-ray period is not the orbital period, but is significantly longer. It is therefore likely that the X-ray period is driven by a super-orbital period, on which, for example, the accretion disk or the neutron star could be precessing. As postulated by Motch et al. (2014) and Hu et al. (2017) it is also possible that the measured X-ray period in fact is a beat period between the orbital period and the super-orbital period, which would imply a super-orbital period on the order of 1500 d. In any case, it seems plausible that the observed 66 d X-ray variability is not a direct tracer of changes in accretion rate, but is also strongly influenced, for example, by the viewing angle onto the accretion disk (similar to NGC 5907 ULX, Fürst et al. 2017).

However, if the X-rays are variable mainly due to geometrical effects and not to a real change in accretion rate, our approach of using the XRT light curve as a tracer for the intrinsic luminosity is flawed. We therefore replace the input light curve with just a simple linear brightening trend, based on the measured data. This approach also results in a very good description of the 2016–2017 pulse period data (Fig. 5), and we find an orbital period of d and an eccentricity of e = 0.05 ± 0.05 (see Table 2). This orbital period is in excellent agreement with the measured optical period.

We also test this model including the 2013 and 2014 period measurements (Fig. 5). However, we cannot achieve an acceptable fit with these epochs. Extrapolating the model back to these times shows deviations up to 300 μs in 2014, but a surprisingly good agreement for the linear trend model in 2013. As we do not have good information about the luminosity between these measurements, different accretion torques might have been acting on the neutron star during the long gaps of no X-ray monitoring. Changes in the luminosity of a factor <2 are enough to change the inferred pulse period by an amount similar to the deviations observed, making this hypothesis likely.

Finally, we compare the measured values of the instantaneous period change Ṗ in each epoch to the values predicted by the model. This information is not used during fitting, and therefore provides an independent test of the model predictions. As shown in Fig. 5d, the model agrees very well with the measured Ṗ values; however, only the NuSTAR measurements are constraining. We show the complete expected Ṗ value, i.e., the combination between secular spin-up and orbital motion, to be directly comparable with the measured values.

Fig. 5. Panel a: Swift/XRT monitoring light curve in the 0.3–10 keV energy band. The three columns focus on epochs 2013, 2014, and 2016–2017. The red curve shows the light curve folded on a 66.9 d period plus a linear brightening trend; the blue dashed line only shows the linear trend. Panel b: measured values of the pulse period as a function of time. Diamonds indicate values taken with NuSTAR and circles values taken with XMM-Newton. The small square shows the period measured in observation 2017A with XMM-Newton, which did not provide a significant detection by itself. The red line is the best-fit model using the pulse profile as input; the blue dashed line a linear trend only. The model takes the secular spin-up as function of intrinsic luminosity and the orbital Doppler motion into account. For details see text. Panel c: residuals as data-minus-model for both models (red: pulse profile input, blue: linear trend only). Panel d: measured and predicted instantaneous change of the period (Ṗ). The XMM-Newton data are shown with gray circles and are not constraining due to the shorter exposure time of the XMM-Newton observations. The models were fitted without taking Ṗ into account. The good agreement between data and model is therefore an independent confirmation of the orbital ephemeris. The time axis is given in MJD.

Open with DEXTER

	Fig. 6. Results of the MCMC run for the orbital determination, based on a linearly brightening input light curve. The blue diamonds give the best-fit values and uncertainties from a standard fit. The contour lines give the 68, 90, and 99% confidence levels, i.e., number of walkers contained within these regions.
Open with DEXTER

3.6. Uncertainty estimation through Monte Carlo simulations

As we are fitting a complex model to a relative sparse dataset, we run Markov chain Monte Carlo (MCMC) simulations to investigate parameter degeneracies and obtain a more realistic estimate of the uncertainties. We use the emcee method implementation in ISIS, which is based upon the parallel “simple stretch” method presented by Foreman-Mackey et al. (2013). We run 200 walkers per free parameter for 5000 steps each. As input light curve we use the simple linear trend, as the model that seems to describe the system’s behavior most realistically. We evaluate the model for the 2016 and 2017 data only.

We analyze the results from the MCMC run after a burn-in period of 30% (1500 steps). Figure 6 shows the results for each parameter pair, together with the best-fit values from a standard χ²-minimization fit. As can be seen in the figure, no strong degeneracies are present, except for those between the argument of periastron ω and the time of periastron passage τ, which can be understood given that the eccentricity is so low. In a circular orbit, these two parameters would be perfectly degenerate within one orbital period. The orbital period is in particular very well determined, and periods longer than 65 d are ruled out at the 99% level.

As degeneracies are small, we also give the estimated 90% uncertainties from the 1D distributions in Table 2. These uncertainties are larger than those estimated from a simple χ² fit, which shows that the main driver of the uncertainties are systematics. In particular the small number of data points we currently have makes a precise determination of all orbital parameters difficult. Nonetheless, we find that the measured orbital period is still significantly different from the X-ray flux period.

Last, we check that there are no acceptable solutions for orbital periods <10 d, which would be more similar to the orbital periods proposed for M82 X-2 and NGC5907 ULX. We run the same analysis, in particular a MCMC simulation with 200 walkers with 5000 steps each in the parameters space 0.5 d ≤ P_orb ≤ 10 d. As expected, we do not find an acceptable solution, with a best fit χ² > 10 and no stable minimum in χ² space. This result is in line with the expectations from the binary properties, where a minimal orbital period of P_orb ≈ 24 d is required for the neutron star to be outside the stellar photosphere for a typical B9Ia star.

4. Discussion and conclusion

We have presented a timing analysis of an XMM-Newton and NuSTAR campaign of the ultra-luminous pulsar NGC7793 P13, carried out in late 2017. Through the variation in the pulse period, we determined an orbital period between ∼63.7 and 64.2 d, consistent with the optical photometric period. In addition to the orbital variation, we observe a strong secular spin-up, which we relate to the transfer of angular momentum from the accreted material. When using a simple linear brightening trend as a function of time and momentum transfer based on the theory by Ghosh & Lamb (1979a), we find a very good description of the pulse period evolution from 2016 to 2017. We are unable to include older data in the model as we lack the necessary information about the intrinsic X-ray flux.

As shown in Fig. 5, we also find a good description of the pulse periods when assuming that the measured X-ray count rate is a direct tracer of the X-ray luminosity. In this case, we find a slightly longer orbital period of 64.2 d. The current available data does not allow us to distinguish between these two scenarios. Another observation campaign in 2019 may be able to clearly distinguish between the models if the average X-ray flux until then is well monitored.

Our best-fit model with a linear brightening trend implies that the observed periodic X-ray variability is not dominated by changes in the accretion rate, but more likely by geometric effects from a precessing accretion disk. This explanation is in line with the fact that the X-ray period is significantly longer than the orbital period, around P_X = 66.9 d and that we observe an almost circular orbit. A possible driver for the disk precession could be the Lense–Thirring effect (Middleton et al. 2018).

It is also possible that the X-ray period is the beat period of the orbital period with a very long super-orbital period. That implies a super-orbital period of around 1500 d. Similar super-orbital periods have been suggested by Motch et al. (2014) and Hu et al. (2017) based on an observed phase jitter of the times of maximum light in the optical. If this interpretation is correct, it implies that the X-ray flux is also significantly influenced by the orbital period. This might be possible through the small but non-negligible eccentricity we find, which could influence the accretion rate as a function of orbital phase. Based on the simple model by Brown & Boyle (1984), variations of a factor of 2 in mass accretion rate over the orbit are certainly possible for an eccentricity of e = 0.1, even for relatively large photospheres of the donor star.

However, the fact that the X-ray period is only slightly longer than the optical period is reminiscent of optical super-hump periods observed in SU Ursae Majoris dwarf novae (Warner 1995) which have also been observed in accreting black-hole systems (O’Donoghue & Charles 1996). The super-hump is caused by a 3:1 resonance between the Keplerian velocity in the outer accretion disk and the orbital period, resulting in tidal forces that produce a finite eccentricity as well as a nodal precession of the accretion disk (e.g., Whitehurst 1988; Lubow 1991). While the P13 system has an inverse mass ratio compared to typical dwarf novae (i.e., the donor is significantly more massive than the accreting object), it is nonetheless extreme (q = M_NS/M_* ≈ 0.08). A similar resonance could therefore take place, observable as periodicity in the X-rays emitted by the accretion disk. Following the calculations by Whitehurst & King (1991), we find that indeed a 4:1 resonance is possible within the tidal radius of the neutron star. To determine the exact super-hump period and strength hydrodynamic simulations are necessary, which are beyond the scope of this work. Additionally, tidal forces due to a misalignment of the neutron star spin axis with the orbital axis could also induce a warp or precession of the accretion disk (Martin et al. 2009).

Using our orbital solution and assuming a canonical mass of about 1.4 M_⊙ for the neutron star and between 18 and 23 M_⊙ for the B9I companion, we can solve the mass function for the inclination. We obtain a semi-major axis of 180–196 R_⊙, assuming P_orb = 64 d and e = 0.1. With the measured a sin i of 208 lt-s (≈90 R_⊙) we find an inclination of around 30°. This is in agreement with the fact that we do not observe eclipses, i.e., we cannot be observing the system edge-on.

With these orbital parameters we find a Roche-lobe radius at periastron of about 103–119 R_⊙. On the other hand, for a temperature of 11 ± 1 kK and an absolute V-band magnitude of M_V = −7.44 we derive a radius of 96–125 R_⊙ for the B9I donor star, with the method used by the genetic algorithm described by Mokiem et al. (2005). The star therefore subtends over about a 60° cone of the sky, as seen from the neutron star. It is very likely that the star fills its Roche lobe even with the small eccentricity implied, allowing for the high accretion rate necessary to power the observed X-ray luminosity of the system. However, we note that the stellar parameters are not very well determined, for example the absolute magnitude could be contaminated by optical emission from the accretion disk. It is therefore possible that the stellar radius is smaller, however, the high observed luminosity strongly argues in favor of Roche-lobe accretion.

P13 shows a hard ULX spectrum (Walton et al. 2018b), consistent with the “hard ultraluminous state” as defined by Sutton et al. (2013). In the “ULX unification model” (Sutton et al. 2013; Middleton et al. 2015) such a hard spectrum implies that we observe the source at relatively low inclination angles, down the evacuated cone surrounded by the super-Eddington disk and its wind (Pinto et al. 2016). If the accretion disk and neutron star inclination is aligned with the orbital inclination, our measurements would therefore be in agreement with this model, assuming the cone has an opening angle of at least 30°. The opening angle (of, at least, the soft X-rays) would need to be larger to also illuminate the star and thereby explain the optical period. However, we note that this model might not hold for accreting neutron stars with strong magnetic fields, and alternative models have been discussed (see, e.g., Koliopanos et al. 2017).

To obtain a significant effect of geometric beaming, which would help to explain the apparent super-Eddington luminosity of P13 without the need to evoke extreme accretion rates, Dauser et al. (2017) find that very narrow cones are necessary. To model the super-orbital variability of NGC5907 ULX (Walton et al. 2016), they require half-opening angles of ≪10°; for half-opening angles of 25° or larger the amplification factor rapidly drops to around 1. We note, however, that the model of Dauser et al. (2017) is based on a rather simplified geometry and does not take the emission geometry of the accretion column into account.

All our orbital solutions require a very small eccentricity e ≤ 0.15. Motch et al. (2014) implied a significantly larger eccentricity between 0.3 and 0.4, based on their modeling of the optical light curve. In particular, they argue that the sharp optical maxima imply an eccentric orbit if they are due to the X-ray illuminated side of the star. However, the accretion geometry of the system is currently unknown, and some of the optical flux could also be contributed by an accretion stream or hot spot on the super-Eddington accretion disk.

For P13, pulsations have so far been seen in every observation that provides a high enough S/N. Two of the other pulsating ULXs have shown pulse drop outs, despite no significant change in flux (Israel et al. 2017a; Bachetti et al. 2014). Here we have shown, however, that P13 also shows strong variability in its pulsed fraction, with values as low as 8% during epoch 2017D. Additionally the measured pulse profile shows considerable variability, deviating from the previously observed sinusoidal shape (Fürst et al. 2016). This variability shows that the accretion and emission geometry changes on relatively short timescales. A detailed analysis of the pulsed spectrum is beyond the scope of this paper, however, and will be presented in a forthcoming work (Walton et al., in prep.).

Acknowledgments

We would like to thank the referee for the very helpful comments that helped to improve the manuscript. We thank M. Kühnel and M. Nowak for the useful discussions. DJW and MJM acknowledge support from STFC Ernest Rutherford fellowships. 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 ASDC (Italy) and Caltech (USA). This research has also made use of data obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States. This work has made use of data supplied by the UK Swift Science Data Centre at the University of Leicester, and also of the XRT Data Analysis Software (XRTDAS) developed under the responsibility of the ASI Science Data Center (ASDC), Italy. This research has made use of a collection of ISIS functions (ISISscripts) provided by ECAP/Remeis observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/).

References

All Tables

All Figures

	Fig. 1. Panel a: light curve from Swift/UVOT, with one data point per observation. Black circles are u-filter, red squares are uvw1-filter, green triangles are uvw2-filter, and blue diamonds are uvm2-filter. Panel b: Swift/XRT light curve in the 0.3–10 keV energy band, with the same binning as the UVOT light curve. Dotted lines mark the times of maximum light in the UVOT data, assuming a constant period of Popt = 63:9 d and dashed lines mark the maximum of the folded XRT profile with a constant period of PX = 66:9 d.
Open with DEXTER
In the text

	Fig. 2. Panel a: epoch-folding (black, left y-axis) and Lomb–Scargle Periodogram (green, right y-axis) of all Swift/UVOT data. Panel b: same as panel a but for the Swift/XRT data. The dotted vertical line marks the strongest period of the UVOT data (which is consistent with the implied orbital period from the pulse period analysis), the dashed line the period of the XRT data. The dashed horizontal line marks the 99.9% false alarm probability of the epoch folding results based on an F-distribution corrected for number of trial periods (Davies 1990).
Open with DEXTER
In the text

	Fig. 3. Example of the epoch-folding in P–Ṗ-space for NuSTAR observation 2017F. The χ2 value at each P–Ṗ pair is color-coded. The graphs at the top and on the right show the cut along the respective dimension for the most significant measurement at P = 415.7050 ms and Ṗ = 6.5 × 10−11s s−1. The orange dashed line shows the 99.99% significance line taking the number of trial periods into account.
Open with DEXTER
In the text

	Fig. 4. Pulse profiles of all available XMM-Newton and NuSTAR observations, folded on the respective measured pulse period. The profiles are repeated once for clarity and phase-aligned by hand to have the minimum at phase 0. For epochs where XMM-Newton and NuSTAR observations where taken simultaneously, NuSTAR is shown in blue using the right-hand y-scale. Each panel is labeled with the epoch name, pulse period, and measured pulsed fraction in percent. For details see text.
Open with DEXTER
In the text

	Fig. 5. Panel a: Swift/XRT monitoring light curve in the 0.3–10 keV energy band. The three columns focus on epochs 2013, 2014, and 2016–2017. The red curve shows the light curve folded on a 66.9 d period plus a linear brightening trend; the blue dashed line only shows the linear trend. Panel b: measured values of the pulse period as a function of time. Diamonds indicate values taken with NuSTAR and circles values taken with XMM-Newton. The small square shows the period measured in observation 2017A with XMM-Newton, which did not provide a significant detection by itself. The red line is the best-fit model using the pulse profile as input; the blue dashed line a linear trend only. The model takes the secular spin-up as function of intrinsic luminosity and the orbital Doppler motion into account. For details see text. Panel c: residuals as data-minus-model for both models (red: pulse profile input, blue: linear trend only). Panel d: measured and predicted instantaneous change of the period (Ṗ). The XMM-Newton data are shown with gray circles and are not constraining due to the shorter exposure time of the XMM-Newton observations. The models were fitted without taking Ṗ into account. The good agreement between data and model is therefore an independent confirmation of the orbital ephemeris. The time axis is given in MJD.
Open with DEXTER
In the text

	Fig. 6. Results of the MCMC run for the orbital determination, based on a linearly brightening input light curve. The blue diamonds give the best-fit values and uncertainties from a standard fit. The contour lines give the 68, 90, and 99% confidence levels, i.e., number of walkers contained within these regions.
Open with DEXTER
In the text