© ESO 2021

1. Introduction

The discovery of pulsations from the ultra-luminous X-ray source (ULX), M 82 X-2 (Bachetti et al. 2014), which led to its identification as an accreting neutron star has opened up a new way of looking at extreme accretion regimes. Such systems, known as ULX pulsars (ULXPs), defy the spherical Eddington limit by orders of magnitude, with the most extreme case being NGC 5907 ULX with luminosities in excess of 10⁴¹ erg s⁻¹ or about 500 times the Eddington luminosity for a standard neutron star (Walton et al. 2016; Israel et al. 2017a; Fürst et al. 2017). One of the most easily studied ULXPs is NGC 7793 P13 (hereafter P13, Fürst et al. 2016; Israel et al. 2017b), as it is nearby (d = 3.40 ± 0.17 Mpc, Zgirski et al. 2017), is isolated from other sources in its host galaxy, and exhibits (almost) persistent pulsations. P13 has a pulse period of around 415 ms and typical luminosities of around 5 × 10³⁹ − 10⁴⁰ erg s⁻¹, clearly placing it in the ULX regime, which is typically defined as L_X > 10³⁹ erg s⁻¹. Fürst et al. (2016) measured a spin-up of Ṗ ≈ −3.5 × 10⁻¹¹ s s⁻¹ which they used to infer a dipole magnetic field of around 1.5 × 10¹² G based on the accretion model of Ghosh & Lamb (1979a).

The source was initially discovered in the X-rays by Read & Pietsch (1999). Motch et al. (2011) identified the companion and mass donor as a B9Ia super-giant. Later, Motch et al. (2014) found 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) interpreted the clearly detected period as the orbital period of the system and find a dynamical mass constraint of < 15 M_⊙ for the compact object, ruling out an intermediate-mass black hole and providing evidence for super-Eddington accretion in the system.

Fürst et al. (2018, hereafter F18) used accurate X-ray period measurements obtained with XMM-Newton and NuSTAR to constrain all parameters of the orbital ephemeris of P13 and found an orbital period of 63.9 ± 0.5 d (statistical uncertainties only), thus confirming the results from Motch et al. (2014). They could also constrain the eccentricity ϵ to be very small (ϵ ≤ 0.14). This almost circular orbit is in slight contradiction with the larger eccentricity implied from the optical light curve, which was necessary to explain the narrow optical maximum under the assumption that the compact object is a black hole (Motch et al. 2014). Updated calculations based on more recent optical data and assuming a neutron star accretor might resolve those differences.

The X-ray flux also shows large variations in a number of different timescales. One important periodic variability is found around 65.05 ± 0.1 d (Hu et al. 2017), based on long-term Swift/XRT monitoring data. This period modulates the flux by a factor of 3−4 during the bright state of P13. Using a longer baseline, F18 updated the results of Hu et al. (2017) and found an X-ray period of 66.8 ± 0.4 d. The difference between the optical and UV as well as the X-ray values might be due resonances in the accretion disk or caused by a warped and precessing accretion disk (Hu et al. 2017, F18). This super-orbital period could also explain the variation in the arrival times of maximum light in the optical (Motch et al. 2014; Hu et al. 2017).

P13 also shows strong long-term X-ray flux variations, for example it exhibits X-ray off-states, where its flux drops below the detection limit of Swift/XRT. On the other hand, it has been in a long, bright X-ray flux state since at least 2016 and likely even since 2013, though we lack dense flux monitoring before 2016. Typical luminosities during this time were around 10⁴⁰ erg s⁻¹. Between January and March 2019 it entered a low state, with the flux dropping drastically over the next few months until it was briefly no longer detectable in individual XRT snapshots (Soria et al. 2020; Hu et al. 2020), before recovering to a low, but significantly detected flux. It is currently unknown if these long-term flux variations are periodic or random.

In July 2020, we obtained a Chandra observation of P13 to obtain a measurement of the low state flux (Walton et al. 2020). We find a luminosity of (4.1 ± 0.5)×10³⁸ erg s⁻¹ (in the 0.3−10 keV band) with a spectrum consistent with the one obtained from the low-state in 2011 and 2012, but with a flux at least an order of magnitude higher¹. Based on contemporaneous Swift/XRT monitoring, it seems that the source had already left the deepest off-state during the Chandra observation. Notably this implies that the current off-state was much shorter than the one in 2011 and 2012, which lasted ≳2 years.

The Chandra data do not provide sufficient time-resolution to measure the pulse period of P13, and are therefore not analyzed here. However, the increased flux encouraged us to ask for further monitoring with XMM-Newton and NuSTAR in July and August 2020. While the flux was higher, neither of those observations yielded detectable pulsations. A full analysis of those data will be presented in a forthcoming publication (Walton et al., in prep.). Swift monitoring through October 2020 shows that P13 continues to be active at a low level.

Here we report on continued NuSTAR and XMM-Newton monitoring, using new data taken in 2018, 2019, and 2020, and following the flux, spectral, and pulse period evolution into the renewed off-state. In Sect. 2 we describe the observations analyzed here and the data reduction methods. In Sect. 3, we discus the data analysis, including the UV and X-ray flux period (Sect. 3.1), the pulse period and its evolution (Sect. 3.2), the evolution of the hardness ratios (Sect. 3.3) and the behavior of the pulsed fraction as a function of time and spectral parameters (Sect. 3.4). We discuss the results in Sect. 4 and conclude the paper with a summary and outlook in Sect. 5. Uncertainties are given at the 90% confidence level, unless otherwise noted.

2. Observations and data reduction

2.1. Swift

Since April 2016, P13 has been monitored by the Neil Gehrels Swift Observatory (Swift; Gehrels et al. 2004) in each of its visibility windows, with a typical cadence of around one week or less and exposure times of 1 ks per snapshot. The visibility constraints result in five observation epochs, each lasting for around nine months (see Fig. 1). Results from previous Swift monitoring data are discussed by F18.

Fig. 1. a: Swift/UVOT light curve in the U-band. b: Swift/XRT lightcurve. The colored lines show the four different assumed flux evolutions, which are used as an input for the pulse period evolution: constant flux (orange), linear brightening and dimming trend (red), measured XRT lightcurve (green), and extrapolated X-ray super-orbital period profile (blue). The dotted line represents the estimated count rate at the Eddington limit. c: pulse period evolution as measured by XMM-Newton (circles) and NuSTAR (diamonds). Superimposed are the best-fit models for the four different input lightcurves in the same colors as in panel b. The brown dotted-dashed lines indicate the times of observations when no pulsations where seen. d: residuals with respect to the linear brightening and dimming input, e: residuals with respect to a constant input, f: residuals with respect to the original lightcurve as input, and g residuals with respect to the X-ray profile input. h: measured (black and gray) and predicted (orange) pulse period derivative Ṗ. The model is based on the constant flux model. i: pulsed fraction in the 3−10 keV energy band. Upper limits are denoted by downward pointing arrows. For details see text.

In addition to the data presented by F18, we extracted 131 XRT (Burrows et al. 2005) observations taken between 2018-04-14 (ObsID 00093149031) and 2020-12-31 (ObsID 00031791109) with the standard Swift/XRT processing pipeline (Evans et al. 2009), thereby extending the data presented by F18 by over three years. The data are binned such that there is a single 0.3−10 keV flux measurement from each observation. Selected observations during the low-state at the end of 2020 were combined manually to yield more stringent upper limits.

We also extracted UVOT (Roming et al. 2005) data from all 131 new observations, following the same method as detailed by F18. In particular, we used a circular source region with a 5″ radius centered on α = 23^h57′50.9″, δ = −32 ° 37′26.6″ and a 15″ circular background region. The data were processed with the corresponding software tasks as distributed by HEASOFT v6.24 and we used uvotsource to extract the source magnitudes.

Figure 1 shows the long-term light curve of these observations obtained with the UVOT (panel a) and XRT (panel b) instruments. These light curves are further discussed in Sect. 3.1.

2.2. NuSTAR

Data from NuSTAR (Harrison et al. 2013) were reduced using the standard pipeline, nupipeline and nuproducts, provided with the NuSTAR Data Analysis Software (v1.8.0), using standard filtering and NuSTAR CALDB v20191219. We extracted source events from both focal plane models (FPMA and FPMB) in circular regions with a radius of 35″ and background events from circular regions with a radius of 120″ on the same detector as the source. We chose the source region size based on optimizing the signal-to-noise ratio (S/N) above 10 keV. Larger regions will add disproportionately more background photons than source photons, reducing the high energy S/N. All time information was transferred to the solar barycenter using the DE-200 solar system ephemeris (Standish et al. 1992). To search for pulsations we combined the source filtered event files for FPMA and B to improve the statistics.

2.3. XMM-Newton

Data from XMM-Newton (Jansen et al. 2001) were reduced with the XMM-Newton Science Analysis System (SAS) v18.0.0, following the standard prescription². We only use data from EPIC-pn (Strüder et al. 2001) in this work, as it provides the necessary fast time resolution to search for pulsations. The data were taken in full frame mode and raw data files were cleaned and calibrated using epchain and transferred to the solar barycenter using the SAS task “barycen” based on the DE-200 solar system ephemeris (Standish et al. 1992).

We extracted source events for all epochs from circular regions with a radius of 40″, following the same method as described in F18. Background spectra were extracted from a source free circular region with a radius of ∼100″, located on the same CCD as P13. We carefully checked all observations for background flaring, but found that it was only problematic for epoch 2017A, which prevents us from measuring the pulse period in that observation (as discussed in F18). See Table 1 for a complete observation log.

3. Analysis

3.1. UV and X-ray periods

Given the much longer timeline of available Swift monitoring data, we updated the long-term periods presented by F18. We used the same approach as presented in F18, that is to say we performed epoch folding (Leahy et al. 1983) and calculated a Lomb-Scargle periodogram (Scargle 1982) for both the Swift/UVOT and the Swift/XRT light curve (Fig. 2). For epoch folding, we used the L-statistics proposed by Davies (1990) for increased sensitivity.

Fig. 2. Results from epoch folding (black, left y-axis) and Lomb-Scargle periodogram (yellow, right y-axis) for the UV light curve (top) and X-ray light curve (bottom). The strongest UV period is marked by the green dotted-dashed line, the strongest X-ray period by the purple dotted line. The new orbital period is indicated by the blue dashed line, the old estimate of the orbital period is shown by the gray dashed line. The shading behind each period indicates their respective uncertainties (using the updated uncertainties for the old orbital period in gray).

Due to the high variability in flux (see Fig. 1b), we neededto normalize the XRT data. F18 used a linear brightening trend and removed it from the data. As such a trend is obviously no longer a good fit, we instead opted to renormalize each epoch to its respective mean count-rate. This approach is the same as used by Hu et al. (2017). No renormalization was done for the UVOT data given their overall stability. Uncertainties (at the 90% level) were determined by simulating 5000 light curves, sampled with the same cadence as the real light curve, with each point drawn randomly from a Poisson distribution based on an interpolation of the respective folded profile.

We find an optical period of d (Fig. 2, top), in very good agreement to previous results (Motch et al. 2014; F18). In the X-rays, we find a period of P_X = 65.31 ± 0.15 d (Fig. 2, bottom), significantly shorter than the 66.8 ± 0.4 d value reported by F18. However, the value we measure here is close to the one presented by Hu et al. (2017): P_X, Hu = 65.05 ± 0.1. Even with this reduction, the X-ray period is very significantly different from the optical period. We checked that the method of removing the underlying variability does not influence the measured value, that is to say we obtain the same results for renormalizing each epoch, subtracting a trend, or not changing the data at all. However, the statistical detection of the X-ray period is significantly improved when using the renormalization for each epoch.

We note that the UV period is much more pronounced during the X-ray low-state in 2020. A continuation of the UV period even during X-ray low-states was already discussed by Motch et al. (2014), who attributed it to the fact that a large precessing accretion disk shields the X-rays from us, but not towards the companion star. The UV variability would then be caused by the X-ray heated side of the companion periodically turning towards us. However, this does not necessarily explain why the UV variability is suppressed during the X-ray high state. A dilution of the UV period due to stronger contribution from the accretion disk to the UV flux seems unlikely, as the average U-band magnitude of the system did not change during the X-ray high state.

While the peak in the X-ray periodogram appears much broader than the peak in the UV periodogram, we do not find any evidence that the X-ray period is quasi-periodic in nature. By splitting the data into smaller parts, we find no indication that the X-ray period is changing in value, however it is more pronounced during the X-ray high state.

3.2. Pulsations

3.2.1. Pulsation search

We searched for pulsations in all new observations (2018A–2020C), using the same accelerated epoch folding search as used in F18. In particular, we searched for pulsations over a grid in the plane defined by the pulse period, P, and its first time derivative, Ṗ, with the data binned into 12 phase bins. To limit the search ranges, we used the secular spin-up and orbital ephemeris found by F18 as an estimator for the expected pulse period during each observation. We then performed a search around that estimated period in a 100 × 100 grid with ΔP = ±0.3 ms and ΔṖ = ±2 × 10⁻⁹ s s⁻¹. Due to their longer duration, the NuSTAR observations provide more constraining measurements (in particular for Ṗ). We therefore performed a second search for the NuSTAR data only, where we zoomed in on the peak found in the previous calculation and searching a 120 × 120 grid with ΔP = ±0.1 ms and ΔṖ = ±6 × 10⁻¹¹ s s⁻¹. We found highly significant pulsations with a significance > 99.5% in all epochs but 2020B and 2020C. Those data were taken during the recovery from the off-state and we discuss them in more detail below. The significance is based on a χ²-statistic corrected for the number of trials corresponding to the bins in the P − Ṗ grid.

To estimate an upper limit on the pulsed fraction in the three epochs where no pulsations are detected (2017A, 2020B, and 2020C), we simulate event files based on an input light curve with the same average count rate as the real data, but with added sinusoidal pulsations, resulting in a pulsed fraction PF_sim. The simulations are based on the method used in the Stingray package (Huppenkothen et al. 2019). For each observation, we simulate event lists with pulsed fractions, PF_sim, 0.1 ≤ PF_sim ≤ 0.9 in steps of 0.01. For each value of PF_sim we generate 100 statistically independent event lists and search each list for pulsations using the epoch folding technique. We do not include Ṗ in these simulations as it does not influence the obtained upper limit. We define the upper limit on the pulsed fraction as the value where at least 90% of all simulated event lists provide a detection of the pulse with at least 99.5% significance. The results are listed in Table 1. We estimated the upper limits assuming a constant pulsed fraction over the whole energy band of the respective instrument.

Uncertainties on P and Ṗ were determined from the extent of the full-width half-maximum (FWHM) contour in the 2D χ² landscape of the epoch folding results. That is, we define the uncertainties in both parameters as the range where the 2D χ² peak has dropped to half of its peak value.

The measured values for P and Ṗ for each observation are given in Table 1 and plotted in Fig. 1c. As can be seen the source continued its secular spin-up with roughly Ṗ ≈ −3.8 × 10⁻¹¹ s s⁻¹ (corresponding to Hz s⁻¹ in frequency space) before pulsations were no longer detected in mid 2020. The measured spin-up in each observation is found to vary on the order of a few ±10⁻¹⁰ s s⁻¹ around the average value, as the orbital Doppler effect dominates there.

3.2.2. Pulse period evolution and orbital ephemeris

We describe the pulse period evolution with a combination of secular spin-up and orbital motion. We apply the same model as described in F18, which allows us to fit for the orbital parameters (orbital period P_orb, eccentricity ϵ, projected semi-major axis a sin i, argument of periastron ω, and time of periastron τ) and requires as input a term related to the accretion of angular momentum. Our first order assumption was that the observed X-ray flux should be an adequate tracer of the accreted angular momentum, following standard accretion theory (Ghosh & Lamb 1979a,b). In this description the pulse period change is expected to be proportional to PL^6/7, where L is the (bolometric) luminosity. As we do not know the exact coupling constant or conversion between observed X-ray count-rates and luminosity, we subsume these conversion in one factor, the spin-up parameter b (for details see Marcu-Cheatham et al. 2015; Bissinger 2016, and F18).

One of the main issues with using the measured X-ray flux as tracer of accreted angular momentum is that our observation history of the X-ray has gaps that can last weeks to months (see Fig. 1b). These are mainly due to gaps in visibility of the source for Swift, and therefore occur roughly once a year. Furthermore the observed X-ray flux may be modulated by intrinsic absorption or changing of the beaming factor of the emitted X-rays, while the actually accreted mass and angular momentum has not changed.

F18 circumvented the problem of the missing data by replacing the measured X-ray flux with two simple models: a linear brightening trend and a variable profile based on folding the Swift/XRT data on the 66.9 d super-orbital X-ray period.

The new data show that a linear brightening trend is no longer a realistic description of the long-term light curve, given the large drop in observed flux in 2019. We therefore modify the trend with a break at around MJD 58300, after which a linear dimming trend is applied (red model in Fig. 1). This approach allows us to build on the solution for the orbit and pulse period evolution found by F18, but also captures the overall shape of the long-term light curve.

However, this model fails to explain the full data-set, leaving large residuals in the X-ray timing data (Fig. 1d). The best fit implies an orbital period of around 65 d (formal uncertainty calculation is not feasible here given the overall bad quality of the fit). Separately, the data before and after January 2018 can be fitted well, however, the best-fit solutions seem to be incompatible with each other, with d and d. Compared to F18, we find slightly larger uncertainties on the orbital period in the data before 2018. Upon closer investigation we found that the uncertainties reported in F18 are underestimated due to a bug in the minimization routine, which has since been fixed.

We also note that the model proposed by Ghosh & Lamb, and in particular the assumptions about how the magnetic field connects to the accretion disk, are likely not applicable in the case of ULXPs. For example, the extreme accretion rates in ULXPs will lead to the formation of geometrically thick accretion disks which were not discussed by Ghosh & Lamb (1979b). Interestingly, Fürst et al. (2016) found that the Ghosh & Lamb theory can explain the spin-up of P13 with a magnetic field of around 1.5 × 10¹² G; however, this only works for the high luminosities observed in 2013−2016. With the lower luminosities observed in 2019, the model predicts much lower maximal spin-up rates, independent of the magnetic field.

A better description of the overall pulse period evolution is obtained when assuming a constant X-ray flux as input (orange model in Fig. 1), that is to say a constant secular spin-up only modulated by the orbital period. This approach implies that the observed X-ray luminosity is not tracing the accretion of angular momentum. This model leaves small residuals around the densely sampled epoch in 2017 (t ≈ 400 d in Fig. 1), however, it provides a much better match to the most recent data during the low flux state of the source. We find an orbital period of around 64.9 d.

We find the same general behavior when comparing a model using the directly measured XRT light curve (green in Fig. 1) as input vs an input based on the super-orbital X-ray profile with a constant average flux (blue in Fig. 1). The large reduction in flux in 2019 in the measured XRT light curve leads directly to an over-prediction of the observed pulse period, while the constant average flux of the profile input provides a much better description of the long-term behavior.

We base our updated orbital calculation on the assumption of a constant spin-up, as it seems to describe the observed observations of the pulse periods best. However, there are still significant outliers in late 2018 (t ≈ 1000 d in Fig. 1) which cannot be explained with this simple model. They are likely caused by brief periods of enhanced accretion, however, they occur at the end of a densely sampled interval, making it unlikely that we missed large X-ray flares that would result in a significant amount of additionally accreted matter and angular momentum. On the other hand, because the X-ray flux is not a good tracer for the amount of accreted angular momentum, it is possible that a spin-up due to enhanced accretion occurred without leaving a measurable trace in the X-ray lightcurve. For calculating the updated ephemeris we therefore first ignore those data points, and discuss the impact of different scenarios to describe them below.

The overall fit of this model is still not very good in terms of χ², with χ² = 64.7 for 7 degrees of freedom (based on 7 orbital parameters and 14 data points). To allow realistic error calculation, which requires a χ² ≈ 1, we add 0.005% of systematic uncertainties on all measurements of the pulse periods (which implies a factor 2−5 increase over the statistical uncertainties and is likely related to timing noise), resulting in a reduced χ², , of 1.06 for the same number of degrees of freedom (d.o.f.). Including the “outlier” data around MJD 58500 results in a best-fit with only a χ² = 88.3 for 8 d.o.f. even with those systematic uncertainties.

Given the complexity of the fit and the low number of degrees of freedom, we also run Markov chain Monte Carlo (MCMC) simulations to estimate the posterior distribution of each parameter. We use an implementation of the “emcee” sampler (Foreman-Mackey et al. 2013) in ISIS, which is based on the method proposed by Goodman & Weare (2010). We use 210 walkers (30 walker per free parameter) and evolve them for 3000 steps. Before calculating the distributions of walkers we use a 20% burn-in period. The results are shown in Fig. 3, together with the best-fit values and uncertainties from the standard χ²-optimizer.

Fig. 3. Orbital parameters distributions based on the results from the MCMC run using a constant spin-up model and ignoring the two measurements around MJD 58500. The blue diamonds and error-bars indicate the results from a standard χ2 optimization. The contours show the 68%, 90%, and 99% confidence intervals.

We find that spin-up strength and initial pulse period (at MJD 57530.0) are very well constrained. We find a best-fit orbital period of d, which is almost a day longer than the orbital period presented by F18 and implied by Motch et al. (2014). The orbital period shows a weak secondary maximum of around 61 d, which also corresponds to a slightly smaller projected semi-major axis and a much larger eccentricity, which seems unphysical and in particular does not describe the densely sampled 2017 data well. We therefore ignore this minimum and report the 1D uncertainties for the orbital period only based on the main peak at 64.87 d.

This longer orbital solution compared to the one presented by F18 is necessary to explain the behavior of the pulse period in late 2019 and early 2020 (t ≈ 1350 d in Fig. 1). These new data strongly constrain the orbital phase, highlighting how important a dense sampling is for constraining the orbital period. With an orbital period of 63.9 d as found by F18, we find that the phase is almost half a period off. While it is possible that the orbital period changes in this system due to loss of angular momentum (see, e.g., Bachetti et al. 2020), the required change would be orders of magnitude larger than expected. We find, however that the older F18 estimate and the updated constraints on the orbital period presented here are still marginally consistent within their ∼99% uncertainties.

The argument of periastron, ω, is basically unconstrained, which is a result of the vanishing eccentricity, ϵ, which is consistent with 0 (similar to the results by F18). Overall, the results from the MCMC run agree well with the values obtained by χ² fitting. We present the 1D uncertainties from the parameter distributions in Table 2.

As mentioned above, this new ephemeris is obtained when ignoring two measurements at the end of 2018. Clearly, the source underwent some stronger spin-up over the course of 2018 than predicted by our model. To test the influence of those data points on our ephemeris, we split the data in two parts, one before January 2018 and one after. We then require that both parts have the same orbital solution, but allow for different spin-up and P(0) values between them. With this, we basically allow a rapid spin-up event at some point during 2018 and possible lower spin-up trend from December 2018 to 2020. We find that the orbital parameters using this model are fully compatible with the values when ignoring the 2018 data. In particular, we find P = 65.05 ± 0.25 d, which is consistent with the orbital period in the previous model and also significantly longer than the UV/optical period.

Regarding the spin-up, we find Ṗ₁ = (−3.93 ± 0.11) × 10⁻¹¹ s s⁻¹ for the first part, and Ṗ₂ = (−3.37 ± 0.13) × 10⁻¹¹ s s⁻¹ for the second part. Both are lower then the best-fit solution presented in Table 2 as this model has an implicit jump of ΔP of around −0.7 ms sometime in 2018. More observations in the future are required to constrain if Ṗ did indeed change in 2018.

3.3. Spectral evolution

In many accreting sources, large changes in flux go together with significant changes in the spectral shape, including X-ray pulsars in the Milky Way (see, e.g., Reig & Nespoli 2013). As many of the spectra show a rather low signal-to-noise ratio (S/N), we restrict ourselves to studying hardness ratios as a proxy for spectral change. A more detailed spectral analysis will be presented in a forthcoming publication (Walton et al., in prep.).

We define three energy bands, soft (S) between 0.5−1.5 keV, medium (M) between 3.0−5.0 keV, and hard (H) between 5.0−10.0 keV. These bands were chosen by eye as they highlight the observed features most clearly, but the exact change of the energy bands does not influence the overall behavior. The soft band is only available for the XMM-Newton data. We measured the flux in each of these bands based on the spectrum in each observation. We define the hardness ratio (HR) as

(1)

where X and Y are the fluxes in the harder and softer energy band, respectively. We plot the hardness ratio as a function of flux in Fig. 4. As can be seen, there is very little variation in the high energy spectrum, with HR(H,M) almost constant over the whole flux range. At lower energies, a slight hardening with increased flux is visible. This could either be due to an increase in absorption or an intrinsic change in the spectral shape (a more in depth analysis of the spectral evolution will be presented in Walton et al., in prep.). We speculate that at higher luminosities, stronger outflows are launched from the super-Eddington accretion disk, which contribute to a larger absorption column.

Fig. 4. Hardness ratio as a function of flux for all XMM-Newton and NuSTAR data. a: HR between the 5.0−10.0 keV and 3.0−5.0 keV energy band for NuSTAR. b: HR between the same bands as in (a) but for XMM-Newton. c: HR between the 3.0−5.0 keV and 0.5−1.5 keV band for XMM-Newton. Data taken in 2020 are shown in red, data from 2019 in blue, and all previous data in black. The dashed lines show the respective average HR for each instrument for all data before 2019. The top y-axis gives the flux as a fraction of the Eddington limit for a 1.4 M⊙ neutron star at a distance of 3.4 Mpc.

During the lowest observed flux (XMM-Newton in epoch 2020C) the source was just around the Eddington limit for a 1.4 M_⊙ neutron star at a distance of 3.4 Mpc (Zgirski et al. 2017). Here the luminosity is estimated based on the 3−10 keV flux, which is roughly a factor of 2 lower than the bolometric (0.5−100 keV) accretion luminosity of P13. At these low flux levels, the top two panels of Fig. 4 suggest a slight softening of the spectrum at higher energies. However, a significant change cannot be claimed given the large measurement uncertainties.

3.4. Pulsed fraction

We calculated the pulsed fraction (PF) in all observations in the 3−10 keV energy band, based on the pulse profile with 12 equally spaced phase bins. We estimate the PF as

(2)

where PP is the pulse profile. The uncertainty of the PF is based on Gaussian error propagation, which is justified as each bin of the pulse profiles contains at least 25 counts.

We find that during the latest NuSTAR observations (epochs 2019B, 2019C, and 2020A) the pulsed fraction was significantly higher than in other observations, reaching up to 60% in the 3−10 keV band. We show the pulsed fraction as function of time in Fig. 1. As shown by F18, the pulsed fraction is typically strongly energy dependent, with higher energies showing higher pulsed fractions. The energy dependence is most significant at low energies (covered XMM-Newton) and levels off at higher energies (covered by NuSTAR). The energy dependence is consistent in most observations, with the exception of the observations in epoch 2019B, which have the highest pulsed fraction overall. In this epoch the pulsed fraction is already very high at low energies and does not show a significant energy dependence.

The pulsed fraction shows an anticorrelation with flux, as shown in Fig. 5, with a Pearson’s correlation coefficient of −0.83 ± 0.07. We estimated the uncertainty on the correlation coefficient via a bootstrapping resampling method using 10 000 iterations. Using the Student’s t-test, we find that the anticorrelation is significant at the > 99.9% level.

Fig. 5. Pulsed fraction in the 3−10 keV band as a function of flux in the same energy band. NuSTAR data are shown as diamonds, XMM-Newton data are shown as circles. Data from 2019 are shown in blue, data from 2020 are shown in red. The brown dotted line shows a simple linear regression to the data. The top y-axis gives the flux as a fraction of the Eddington ratio for a 1.4 M⊙ neutron star at a distance of 3.4 Mpc.

On the other hand, we do not find a strong correlation overall between the spectral shape, as measured by the hardness and the pulsed fraction (Fig. 6). When taking only the 2019 and 2020 data into account, a correlation can be implied, though it is not statistically significant.

Fig. 6. Pulsed fraction in the 3−10 keV band as a function of the hardness ratio between the 5.0−10.0 keV and the 3.0−5.0 keV energy band. NuSTAR data are shown as diamonds, XMM-Newton data are shown as circles. Data from 2019 are shown in blue, data from 2020 in red.

4. Discussion

We have presented an analysis of the X-ray pulsations seen from NGC 7793 P13 between 2016 and 2020. During this period, the source was mostly in an active state, and showed a constant long-term spin-up. However, in 2019 the observed flux faded significantly, dropping below the detection threshold for our Swift/XRT monitoring, before recovering to a more stable (but still low) flux level. Despite this flux evolution, our X-ray timing results imply that the long-term spin-up continued at a similar rate to that seen in the high-flux state.

In addition to tracking the timing data, we have also explored whether the strength of the pulsed signal evolves with both flux and spectral hardness of P13. The pulsations appear to be strongest at low fluxes, but we find little evidence for any dependence on spectral shape. These results also allow us to make a preliminary assessment of the spectral evolution of P13 across this period. Although during the high-flux period we find little evidence for large spectral variations, we do see some interesting hysteresis associated with the more recent flux evolution.

Using a constant spin-up approximation we have updated the orbital ephemeris of NGC 7793 P13 and find values inconsistent with the ones presented by F18. In particular, we find an orbital period of 64.86 ± 0.19 d based on extensive MCMC simulations. This period is larger than the best-fit orbital period presented by F18, and also longer than the periodicity seen in the UV. On the other hand, it is very close to the revised X-ray period, which we find is 65.31 ± 0.15 d. It is currently unclear how to interpret these different periods in a physical context, and it is particularly puzzling how the optical flux seems to vary on time-scales faster than the orbital period.

It is possible that our estimate of the orbital period has larger systematic uncertainties than implied. As discussed, not all measured periods fit the curve well; in particular, the two measurements around MJD 58450 cannot be reconciled with any simple spin-up model. Hence it is possible that there is timing noise present or that there are unobserved spin-up or -down episodes that we cannot model. Adding ad-hoc flares in the gaps of the XRT monitoring, it is possible to find a model describing these points well, under the assumption that this modified X-ray flux is related to the spin-up value. However, we still find that an orbital period of around 65 d is required to describe all the data and this ad-hoc flux evolution is not based on any observational evidence.

The decoupling between the observed spin-up and the X-ray flux could indicate that strong obscuration occurs during the drop in flux in 2019, while the intrinsic accretion is continuing unabated (or shows flares that result in short-term spin changes). This behavior and scenario is similar to the one proposed for NGC 300 ULX-1, another highly variable ULXP (Vasilopoulos et al. 2019). However, if absorption and obscuration is the reason for the diminishing X-ray flux, we would expect to see a significant hardening of the observed X-ray spectrum, which is not the case (Fig. 4). In fact, we find rather the opposite behavior: the source is getting softer at lower fluxes.

We find a clear anticorrelation between the pulsed fraction and the source flux in the 0.5−20 keV energy band, with pulsed fractions as high as 60% during the low states in 2019 and early 2020. This anticorrelation is clearly present even outside the lowest fluxes. It seems to indicate that at lower fluxes the accretion column, which is responsible for the pulsed flux, dominates. According to Walton et al. (2018), the pulsed flux of the accretion column can be described by a power law with an exponential cut-off at high energies. Typically this component dominates at higher energies, with non-pulsed emission likely associated with the accretion disk also seen at lower energies. These results may suggest a lower relative contribution in the observed bands from the disk. This will be explored in more detail in future work.

On the other hand, the pulsed fraction does not show a significant correlation with spectral hardness (Fig. 6). We would expect a strong correlation if indeed the pulsed fraction increases because the hard accretion column starts to dominate. Instead, the change pulsed fraction might be related to a changed scattering time within the cone of the accretion disk and wind. This cone confines the emitted X-rays to its opening angle, causing so-called “beaming”. Before photons emitted from the neutron star can escape this cone, they might undergo a number of scatterings, causing a significant delay in their arrival time. For a large enough cone this might lead to a smeared out pulse profile with a lower pulsed fraction. We expect large accretion disk cones and larger beaming fractions at higher luminosities, providing a possible way to explain the correlation between pulsed fraction and flux without a significant change in the spectral shape.

5. Conclusion and outlook

NGC 7793 P13 continues to surprise us with new behavior. It is one of only two ULXPs for which the companion star is identified (the other one being NGC 300 ULX-1, Heida et al. 2019, while NGC 1313 X-2 has a known optical counterpart, but the origin of the optical emission is not yet identified, Sathyaprakash et al. 2019; Grisé et al. 2008), and the only ULXP for which the full ephemeris can be determined. However, the details of this ephemeris are still unclear. With the most recent data the best-fit orbital period is 64.9 d almost a day longer than the optical and UV period, and about 0.5 d shorter than the X-ray period. Further observations of the pulse period evolution will allow us to obtain a better understanding if this difference in periods is real or due to a systematic effect in our measurement.

We have also found a correlation between the flux and the pulsed fraction and have shown that the pulsed fraction can change significantly without a measurable change in spectral shape. Forthcoming detailed spectral modeling (Walton et al., in prep.) will allow us to investigate this behavior in more detail and probe different scenarios of obscuration by neutral or highly ionized material. In addition, continued measurement of the pulsed fraction at different flux levels will allow us to fill in the parameter space and investigate if clear changes in accretion geometry occur at certain fluxes.

A major step forward in our understanding of ULXPs would be provided by updated models no how the torque of the accreted material is transferred onto the neutron star and how the magnetic field couples with the accretion disk in the case of geometrically thick, super-Eddington accretion disks.

Acknowledgments

We would like the thank the referee for the very useful comments that helped to improve the manuscript. 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 the 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 made use of data supplied by the UK Swift Science Data Centre at the University of Leicester, and also made use 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/). The material is based upon work supported by NASA under award number 80GSFC17M0002.

References

All Tables

All Figures

Fig. 1. a: Swift/UVOT light curve in the U-band. b: Swift/XRT lightcurve. The colored lines show the four different assumed flux evolutions, which are used as an input for the pulse period evolution: constant flux (orange), linear brightening and dimming trend (red), measured XRT lightcurve (green), and extrapolated X-ray super-orbital period profile (blue). The dotted line represents the estimated count rate at the Eddington limit. c: pulse period evolution as measured by XMM-Newton (circles) and NuSTAR (diamonds). Superimposed are the best-fit models for the four different input lightcurves in the same colors as in panel b. The brown dotted-dashed lines indicate the times of observations when no pulsations where seen. d: residuals with respect to the linear brightening and dimming input, e: residuals with respect to a constant input, f: residuals with respect to the original lightcurve as input, and g residuals with respect to the X-ray profile input. h: measured (black and gray) and predicted (orange) pulse period derivative Ṗ. The model is based on the constant flux model. i: pulsed fraction in the 3−10 keV energy band. Upper limits are denoted by downward pointing arrows. For details see text.

In the text

Fig. 2. Results from epoch folding (black, left y-axis) and Lomb-Scargle periodogram (yellow, right y-axis) for the UV light curve (top) and X-ray light curve (bottom). The strongest UV period is marked by the green dotted-dashed line, the strongest X-ray period by the purple dotted line. The new orbital period is indicated by the blue dashed line, the old estimate of the orbital period is shown by the gray dashed line. The shading behind each period indicates their respective uncertainties (using the updated uncertainties for the old orbital period in gray).

In the text

	Fig. 3. Orbital parameters distributions based on the results from the MCMC run using a constant spin-up model and ignoring the two measurements around MJD 58500. The blue diamonds and error-bars indicate the results from a standard χ2 optimization. The contours show the 68%, 90%, and 99% confidence intervals.
In the text

Fig. 4. Hardness ratio as a function of flux for all XMM-Newton and NuSTAR data. a: HR between the 5.0−10.0 keV and 3.0−5.0 keV energy band for NuSTAR. b: HR between the same bands as in (a) but for XMM-Newton. c: HR between the 3.0−5.0 keV and 0.5−1.5 keV band for XMM-Newton. Data taken in 2020 are shown in red, data from 2019 in blue, and all previous data in black. The dashed lines show the respective average HR for each instrument for all data before 2019. The top y-axis gives the flux as a fraction of the Eddington limit for a 1.4 M⊙ neutron star at a distance of 3.4 Mpc.

In the text

Fig. 5. Pulsed fraction in the 3−10 keV band as a function of flux in the same energy band. NuSTAR data are shown as diamonds, XMM-Newton data are shown as circles. Data from 2019 are shown in blue, data from 2020 are shown in red. The brown dotted line shows a simple linear regression to the data. The top y-axis gives the flux as a fraction of the Eddington ratio for a 1.4 M⊙ neutron star at a distance of 3.4 Mpc.

In the text

	Fig. 6. Pulsed fraction in the 3−10 keV band as a function of the hardness ratio between the 5.0−10.0 keV and the 3.0−5.0 keV energy band. NuSTAR data are shown as diamonds, XMM-Newton data are shown as circles. Data from 2019 are shown in blue, data from 2020 in red.
In the text