M51 ULX7 is among a small group of known ultraluminous X-ray pulsars (ULXPs). The neutron star powering the source has a spin period of 2.8 s, orbits its companion star with a period of 2 days, and a superorbital period of 38 days is evident in its X-ray lightcurve. Here we present NuSTAR and XMM-Newton data on the source from 2019 obtained when the source was near its peak brightness. We detect the pulsations, having spun up at a rate of 3 ± 0.5 × 10−10 s s−1 since they were previously detected in 2018. The data also provide the first high-quality broadband spectrum of the source. We find it to be very similar to that of other ULXPs, with two disk-like components, and a high-energy tail. When combined with XMM-Newton data obtained in 2018, we explore the evolution of the spectral components with superorbital phase, finding that the luminosity of the hotter component drives the superorbital flux modulation. The inclination the disk components appear to change with phase, which may support the idea that these superorbital periods are caused by disk precession. We also reexamine the superorbital period with 3 yr of Swift/XRT monitoring, finding that the period is variable, increasing from 38.2 ± 0.5 days in 2018–2019 to 44.2 ± 0.9 days in 2020–2021, which rules out alternative explanations for the superorbital period.
1.Introduction
M51 ULX7 was first detected as an X-ray source in the galaxies of M51 (NGC 5194/5) by the Einstein X-ray Observatory (Palumbo et al. 1985). The source was observed with an X-ray luminosity of ∼1039 erg s−1, which the authors noted was brighter than any X-ray source in our own Galaxy, or M31, and that it exceeded the Eddington luminosity of a 1 M⊙ object. They suggested this X-ray source, and others like it, could be powered by a neutron star experiencing super-Eddington accretion in a nonspherically symmetric accretion flow geometry or could indicate the presence of a more massive black hole. These luminous X-ray sources would then come to be known as ultraluminous X-ray sources (ULXs; Fabbiano 1989; see also recent reviews by Kaaret et al. 2017 and Fabrika et al. 2021). Thanks to the sensitivity and spatial resolution of XMM-Newton and Chandra, we now know hundreds of ULXs (e.g., Liu & Mirabel 2005; Swartz et al. 2011; Walton et al. 2011; Earnshaw et al. 2019; Kovlakas et al. 2020).
While the assumption that ULXs were powered by black holes gained the most traction since their discovery by Einstein, the neutron star hypothesis was eventually proven, at least for a handful of sources. This began in 2014 with the detection of coherent pulsations from the ULX M82 X-2 by the Nuclear Spectroscopic Telescope Array (NuSTAR) observatory (Bachetti et al. 2014), determining the source to be powered by a neutron star, since black holes are incapable of producing such signals. This was followed by NGC 5907 ULX (Israel et al. 2017a), NGC 7793 P13 (Israel et al. 2017b; Fürst et al. 2017), NGC 300 ULX (Carpano et al. 2018), NGC 1313 X-2 (Sathyaprakash et al. 2019), and the subject of this paper, M51 ULX7 (Rodríguez Castillo et al. 2020). SMC X-3 (Tsygankov et al. 2017), Swift J0243.6 + 6124 (Wilson-Hodge et al. 2018) and RX J0209.6-7427 (Chandra et al. 2020; Vasilopoulos et al. 2020b) also briefly became ULX pulsars, and another candidate pulsating ULX was also recently reported in NGC 7793 (Quintin et al.2021).
From timing analysis of XMM-Newton data, Rodríguez Castillo et al. (2020) determined the spin period of M51 ULX7 to be 2.8 s, and that the neutron star was in a 2 day orbit with a >8 M⊙ companion star, making it a high-mass X-ray binary system (HMXB). The long-term, secular spin-up rate was also found to be ∼10−9 s s−1. Hu et al. (2021) and Vasilopoulos et al. (2021) also found evidence for periodic dips in the Chandra X-ray lightcurve that are associated with the 2 day binary orbital period which they interpret as eclipses, implying that the orbit of the neutron star and its donor star is seen at high inclination.
The spin period is similar to the ∼1 s spin periods of M82 X-2, NGC 5907 ULX, and NGC 7793 P13. All four sources also share the common characteristic of having periodic flux modulations in their long-term X-ray lightcurves, with the latter over a range of 60–80 days (Walton et al. 2016; Fürst et al. 2018; Brightman et al. 2019). The flux modulation from M51 ULX7 was found to have a period of 38 days from Swift/XRT monitoring (Vasilopoulos et al. 2020a; Brightman et al. 2020b). These periodic flux modulations, which in most cases are longer than the orbital period of the system, have been interpreted as precession of a large-scale height disk, possibly caused by the Lense–Thirring effect (e.g., Middleton et al. 2018). However, the origin of these superorbital periods is still a matter of debate.
While the aforementioned ULXs have now been determined to be powered by a neutron star, it is still unknown what fraction, if any, are powered by black holes. From an X-ray spectral standpoint, the neutron-star-powered ULXs appear very similar to ones with unknown accretors, albeit among the hardest, implying that the vast majority could be powered by neutron stars (Koliopanos et al. 2017; Pintore et al. 2017; Walton et al. 2018a; Gúrpide et al. 2021). The spectral shape consists of two disk-like components, a cooler one which may come from the outer regions of an accretion disk or the photosphere of an outflow (Qiu & Feng 2021), and a hotter component, which may originate from the inner regions of the accretion disk (Walton et al. 2018a), an accretion curtain (Mushtukov et al. 2017), or Compton up-scattering (Titarchuk 1994). A high-energy tail is also seen when NuSTAR data are available, and appears to be associated with the pulsed component (Walton et al. 2018a).
While NuSTAR has observed M51 ULX7 on two previous occasions with lower energy coverage, the exposure time was either too short for a good quality spectrum (Earnshaw et al. 2016), or the source was caught at a low flux (Brightman et al. 2018). In NuSTAR Cycle 5, we obtained joint NuSTAR and XMM-Newton observations of M51 ULX7, timed to occur at the peak of the periodic flux modulation, with the aims of obtaining a high-quality broadband X-ray spectrum of ULX7, modeling its emission components, and tracking its pulsations. We present the results from these observations in this paper. We assume a distance of 8.58 ± 0.10 Mpc to M51, derived from the tip of the red giant branch method (McQuinn et al. 2016).
2.X-Ray Data Reduction
2.1.NuSTAR
NuSTAR (Harrison et al. 2013) observed M51 from 2019 July 10 (UT)05:56:09 to 2019 July 14 (UT)01:26:09 with an exposure of 169 ks (obsID 60501023002). We used heasoft v6.28, nustardas v2.0.0 and caldb v20201101 to analyze the data. We produced cleaned and calibrated events files using nupipeline with the settings saacalc=3 saamode=OPTIMIZED tentacle=yes to account for enhanced background during passages of the South Atlantic Anomaly (SAA), which reduced the exposure time to 162 ks. We used nuproducts to produce spectral data, including source and background spectra, and response files. A circular region with a radius of 30'' was used to extract the source spectra. A circular region with a radius of 100'' was used to extract the background spectra, taking care to extract the background from the same chip as the source. For timing analyses, we used the heasoft tool barycorr to apply a barycentric correction to the event times of arrival.
2.2.XMM-Newton
XMM-Newton (Jansen et al. 2001) observed M51 from 2019 July 11 (UT)10:47:26 to 2019 July 12 (UT)08:09:24 with an exposure of 77 ks (obsID 0852030101). We used xmmsas v18.0.0 to analyze the data. We first identify periods of high background by creating a lightcurve of the events in the 10–12 keV band, creating good time intervals where the rate was less than 0.1 counts s−1 in this band, leaving 69 ks of data. Events were selected with PATTERN ≤ 4 for the pn and PATTERN ≤ 12 for the metal oxide semi-conductor (MOS) cameras. A circular region with a radius of 30'' was used to extract the source spectrum. A circular region with a radius of 60'' was used to extract the background spectra, on the same chip as the source and also in the galaxy in order to account for the soft diffuse emission the source is embedded in (e.g., Earnshaw et al. 2016). Data from the pn and both MOS instruments were extracted in this way. For timing analyses, we used the xmmsas tool barycen to apply a barycentric correction to the event times of arrival.
2.3.Swift
We used the online tool provided by the University of Leicester 8 (Evans et al. 2007, 2009) to extract the Swift/XRT (Burrows et al. 2005) lightcurve of ULX7. All products from this tool are fully calibrated and corrected for effects such as pile-up and the bad columns on the CCD. We selected observations with target IDs 11417, 30083 and 32017, and binned the lightcurve in time, with a maximum bin size of 500 ks (5.79 days) and a minimum detection of 2.5σ. The lightcurve is plotted in Figure 1, with the time of the 2019 NuSTAR and XMM-Newton observation marked, as well as the 2018 XMM-Newton observations.
3.Pulsation Analysis
We searched for pulsations using the fast Z2 n search implemented in the hendrics tool HENzsearch (Bachetti 2015). This tool folds the data along a grid of frequencies and frequency derivatives and calculates the Z2 n statistic starting from the folded profiles (see Huppenkothen et al. 2019; Bachetti et al. 2021). We used n = 1 (Rayleigh test) and folded profiles of 16 bins (adequate for using the 
search with the binned approximation).
We found a clear peak in the f–
plane in the XMM-Newton data; however, no such signal was found in the lower-count NuSTAR data. We noted the best solution from HENzsearch and used HENphaseogram to calculate the pulse time-of-arrival (TOA) in 16 intervals during the observation. Following this, the graphical tool pintk in pint (Luo et al. 2021) was used to fit the TOAs and get a timing solution and its uncertainty. We find f = 0.358842(5) Hz (P = 2.78674 s), and 
. We used the solution from pint to phase tag the events using HENphasetag.
We did not find evidence for the orbit of the neutron star and its companion when searching for a second derivative in the pulse frequency. This is likely due to too few counts and only a fraction of the orbit being covered by the XMM-Newton observation. Therefore we could not make a correction for the orbit in the determination of the spin period. This means that the observed first derivative in the pulse frequency is the sum of both the secular spin up of the neutron star and the orbital motion. Rodríguez Castillo et al. (2020) calculated the maximum delay/advance introduced when not correcting for the orbital motion, which was of the order of 1 ms. In the following, when inferring the secular spin up of the pulsar, we assume this 1 ms value to be the absolute uncertainty of the period.
Using the above timing solution and the phase-tagged events files, we create pulse profiles by binning the lightcurve in 16 equally sized phase bins, and energy bins of 0.5–0.1 keV, 1–2 keV, and 2–10 keV. We plot these in Figure 2. Using these, we also calculate the pulse fraction in each bin, defined as the amplitude of the pulse divided by the mean count rate. These are plotted in Figure 3. Also shown are the pulse fractions as determined by Rodríguez Castillo et al. (2020) from the 2018 XMM-Newton observing campaign. We derive an upper limit of 74% to the pulse fraction in the 10–20 keV band from the NuSTAR data by searching 170 Fourier frequencies, and assuming 
.
4.Spectral Analysis
We begin analyzing the 2019 NuSTAR Cycle 5 and XMM-Newton spectral data on ULX7 by grouping the spectra with a minimum of 1 count per bin. We load the FPMA, FPMB, pn, MOS1 and MOS2 source spectra into xspec, subtracting the background, and consider energies 3–20 keV for the NuSTAR data (the background dominates the source above 20 keV so we do not consider these data), and 0.2–10 keV for the XMM-Newton data. We use the C-statistic, suitable for the low number of counts per bin here, and we use a constant term to take account of cross-calibration uncertainties between instruments, which are typical of those found in Madsen et al.(2015).
Rodríguez Castillo et al. (2020) carried out fits to the series of high-quality XMM-Newton spectra they obtained in 2018 when the pulsations were detected. Their best-fit model consisted of an absorbed disk blackbody model plus a hotter blackbody model (tbabs∗ztbabs∗(diskbb+bbodyrad)), which we fit here. The tbabs model accounts for absorption in our Galaxy, fixed at 3.19 × 1020 cm−2 (HI4PI Collaboration et al. 2016) and the ztbabs accounts for absorption at the redshift of M51, z = 0.002, left as a free parameter. This model resulted in C = 4027.81 with 4286 d.o.f.
However, there appears to be an excess at energies above 10 keV when fitting with this model to the 2019 NuSTAR+XMM-Newton data (Figure 4). Rodríguez Castillo et al. (2020) noted that in Walton et al. (2018a) a third component, a power-law with a high-energy cut off (cutoffpl), is used to model this component, which is attributable to the accretion column. This pulsed component was isolated in the ULX pulsars M82 X-2, NGC 7793 P13 and NGC 5907 ULX by extracting spectra from the brightest and the faintest quarters (Δϕpulse = 0.25) of the pulse cycle and subtracting the latter from the former (i.e., pulse on–pulse off; Brightman et al. 2016; Walton et al. 2018a, 2018b).
We do the same pulse on–pulse off analysis here for the XMM-Newton data on ULX7 using the timing solution of the pulses found in Section 3, and fit it with tbabs∗ztbabs∗cutoffpl, where the absorption components have been fixed to the values from the time-averaged spectrum. We cannot constrain Ecut due to the lack of high-energy coverage, so we fix Ecut = 8.1, the average values from Walton et al. (2018a). We can constrain Γ = 0.8 ± 0.3 with C = 881.89 with 946 d.o.f., which is consistent with the average value of 0.5 from Walton et al. (2018a). We plot this spectrum in Figure 5.
We then add this component to the time-averaged spectral fit, keeping Γ and Ecut fixed, but leaving the normalization free. This results in C = 3954.59 with 4285 d.o.f., presenting an improvement of the overall time-averaged fit.
We note that if a diskpbb model were used for the hottest disk component instead of bbodyrad, consistent with the models used in Walton et al. (2018a), we get a better fit of C = 3951.17 with 4284 d.o.f. Therefore, we present this diskbb+diskpbb+cutoffpl model as our best fit. We present the results from the spectral modeling described above in Table 1.
5.Evolution of the Spectral Parameters
With the newly determined best-fit model found in Section 4, we proceed to apply it to the high-quality data from XMM-Newton in 2018 in order to explore the evolution of the X-ray spectral parameters with luminosity and superorbital phase. While these data were presented by Rodríguez Castillo et al. (2020), they used a slightly different model, so we cannot strictly compare their results to ours. We note that these models are purely phenomenological, and any physical interpretations should have this added caveat. We present the details of all the NuSTAR and XMM-Newton observations used in this work in Table 2.
We reduce the XMM-Newton obsIDs 0824450901, 0830191401, 0830191501, and 0830191601 in the same way as described in Section 2.2. We then fit the spectra with the diskbb+diskpbb+cutoffpl model described above. As found by Rodríguez Castillo et al. (2020), ULX7 was seen in a low-flux state in obsID 0830191401, and only the low-temperature diskbb can be seen. We remove the cutoffpl model from the fit since it clearly over fits the data, and we fix the temperature and p parameter of the diskpbb model to typical values in order to place an upper limit on the normalization and flux of this component.
For each epoch, we calculate the total intrinsic luminosity and the luminosity of the two disk components separately, determined using the model component cflux and calculated in the 0.3–10 keV range. To get the superorbital phase, we divide the number of days since 2018 May 1 by the period of 38.3 days, and subtract off the nearest integer. We also adjust the phase so that phase = 0 corresponds to the peak of the superorbital modulation. In Brightman et al. (2020b), we showed that the superorbital period of 38.3 days was consistent across the 2018–2019 epoch we study here.
First we plot the total intrinsic luminosity, and the luminosity of the two disk components as a function of superorbital phase, and compare these to the average total luminosity as seen from the Swift/XRT monitoring. Figure 6 shows that the total luminosity measured by XMM-Newton agrees in general with that of the average luminosity seen by Swift/XRT. In terms of the two disk components, the cool component's luminosity remains relatively constant, whereas the hot component's luminosity varies more, albeit with large uncertainties.
We also plot the pulsed fraction derived in Section 3 and in Rodríguez Castillo et al. (2020), and how it varies with superorbital phase, in Figure 6. We see a similar picture here, where the pulsed fraction in the low-energy band varies relatively little over the superorbital period, but the higher-energy pulsed fraction varies strongly. This is likely driven by the variation in the nonpulsed, hot disk component.
We then proceed to plot the spectral parameters of the disk components against the intrinsic luminosity and superorbital phase in Figure 7. We do not see any statistically significant correlations between the spectral parameters and LX for either the cool diskbb component or hot diskbb component. We do see that NH, Tin, cool, Ncool, and Nhot appear to vary with the phase of the superorbital period, however we likely do not have enough data to claim a significant relationship. If the dependence is real, then since the normalization of the disk components depends on the orientation of the disk, this may imply that the inclination of the disks is changing with superorbital phase, or, in other words, the disks are precessing.
However, since the disk temperature and normalization parameters are degenerate with each other, it is possible that these variations are driven by the degeneracy instead. We investigate this possibility by using Markov Chain Monte Carlo methods to map out the Tin, cool and Ncool parameter space. In xspec we use the Goodman–Weare algorithm with eight walkers and a total length of 10,000 steps with a burn-in phase of 5000 steps. Figure 8 shows the results of this. While a degeneracy can be seen between the two parameters, the results map out regions which are almost mutually exclusive, implying that the degeneracy is not the cause of the possible dependence of these parameters on superorbital phase.
6.Evolution of the Superorbital Period
In Brightman et al. (2020b), we showed that ULX7 exhibited a superorbital period of 38 days over 500 days (∼1.5 yr) from 2018 May (see also Vasilopoulos et al. 2020a). However, when we applied our analysis to the most recent, full ∼3 yr data set, the signal appears much weaker, with a peak in the L statistic at 34.1 days and a peak in the Lomb–Scargle periodogram at 39.0 days. This implies that after 500 days the superorbital period disappears, or changes phase and/or period. Figure 9 shows the full 3 yr lightcurve with the average profile of the flux modulations from the first 500 days overplotted. The deviation from the profile after 500 days is clear in the residuals, and exhibits both positive and negative deviations and therefore cannot be explained by the anomalous low-flux states seen by Vasilopoulos et al. (2020a).
We investigate further by splitting the lightcurve into smaller sections in time to determine how the superorbital period evolves. As in Brightman et al. (2020b), we use epoch folding on the L statistic (Davies 1990), using 10 phase bins. This time, however, since we know the approximate period, we search over a narrower range of 30–50 days in 200 equally spaced bins for the epoch folding. We find that a 300 day section is sufficient to recover the 38 day period at the beginning of the lightcurve, and then we progress through the lightcurve with steps of 30 days.
In Brightman et al. (2020b), we used simulations to determine the false alarm rate for the 38 day period signal, finding it to be >99.9% significant. We repeat this analysis for our 300 day bins, simulating 10,000 lightcurves with 2000 s resolution and a red noise power spectrum, and we sample them with the same observational sampling as the real lightcurves. We then note the largest peak in each periodogram, irrespective of period. We define the false alarm rate as the number of simulated lightcurves that produce a peak as high as the real one, divided by the total number of simulations. We do this for each of the time bins since each time bin has a different number and spacing of observations. For each time bin, we determine the L-stat level, which corresponds to a 0.3% false alarm rate, equivalent to a 3σ detection. In Figure 10 we plot the superorbital periods that are detected at >3σ against time.
In order to determine the uncertainties on the periods detected above, we use a Monte Carlo technique to resample the lightcurve, where the count rate of each observation is a random number that is taken from a Gaussian distribution with a mean equal to the observed count rate and a standard deviation equal to the uncertainty on the count rate. We do this 10,000 times, rerunning our analysis on each lightcurve, and noting the peak period. We define the 1σ uncertainty as the standard deviation of periods recovered. We plot these derived uncertainties in Figure 10.
We find that at the beginning of the lightcurve the superorbital period is detected at >3σ at ∼38 days, as found in Brightman et al. (2020b), but at the end of the lightcurve it is detected at ∼44 days. The highest peaks in the periodograms are at 0–300 days, with P = 38.2 ± 0.5 days, and 780–1080 days where P = 44.2 ± 0.9 days. These time bins are independent and the findings indicate that the superorbital period has increased in length. There are no obvious correlations between the superorbital period and the number of observations per bin, or the average count rate.
While the above approach accounts for the statistical uncertainties on the individual measurements, it does not account for possible effects due to sampling. In order to do this we employ a bootstrapping method, whereby we randomly exclude 10% of the observations in each time bin. Repeating this 10,000 times, we find the standard deviations on the periods derived were 0.13 for the 38.2 day period 0.58 for the 44.2 day period. If we increase the number of excluded pointings to 20%, this becomes 0.28 and 1.83, respectively. This implies the statistical uncertainty on the measurements dominates the uncertainty in the superorbital period.
Finally, we also check if the longer superorbital period seen at 780–1080 days can be produced by the 38.2 day signal, but with the lightcurve sampling at 780–1080 days. We do this by projecting the 38.2 day profile to 780–1080 days and simulating a lightcurve with the same observing strategy during that time; again, we do this 10,000 times. However, none of these simulations produce a peak at 44.2 days or greater with a height at least equal to the observed periodogram. It is therefore highly unlikely that the Swift observing strategy caused the observed change in superorbital period.
We show the lightcurves and folded lightcurves from the 0–300 days and 780–1080 days epochs in Figure 11. The profiles from the two lightcurves appear similar, albeit perhaps with the most recent one being more peaked and less sinusoidal.
7.Discussion
7.1.The Spin of the Neutron Star
The pulsations detected with P = 2.78674 s in our 2019 XMM-Newton data set of M51 ULX7 represent a secular spin up of the neutron star from P = 2.79812 s in 2018 (Rodríguez Castillo et al. 2020), with an average spin-up rate of 3 × 10−10 s s−1 over the 425 days spanning the observations. As described in Section 3, we were not able to correct the 2019 data for the orbital motion of the pulsar. We determined that this produced an uncertainty of order 1 ms in the period. This then leads to an uncertainty of 0.5 × 10−10 s s−1 in the derived secular spin-up rate. The secular spin-up rate is similar to that calculated for the short timescale covered by the 2018 data, and appears comparable with the spin up of ∼10−9 s s−1 calculated from observations spanning 13 yr, which is in the range of other ULX pulsars (Rodríguez Castillo et al. 2020).
7.2.The Broadband X-Ray Spectrum
We presented here the first high-quality broadband X-ray spectrum with coverage above 10 keV of M51 ULX7, obtained with NuSTAR and XMM-Newton, over the energy range 0.2–20 keV. Other members of the known ULX pulsar population with at least one high-quality broadband X-ray spectrum are NGC 5907 ULX (Walton et al. 2015), NGC 7793 P13 (Walton et al. 2018b), and NGC 300 ULX (Carpano et al. 2018). M82 X-2 has broadband spectra but with significant confusion due to the proximity of M82 X-1 (Brightman et al. 2020a), and NGC 1313 X-2 is not well detected above 10 keV (Bachetti et al. 2013). Therefore adding to this small sample of ULX pulsars with a high-quality broadband X-ray spectrum is significant. Indeed, the sample of ULXs with high-quality broadband spectra in general, be it with a known neutron-star accretor or an unknown accretor, is also small (Walton et al. 2018b).
The broadband X-ray spectrum of M51 ULX7 can be well described by two disk-like components, plus a higher-energy component associated with the pulsations. This is qualitatively very similar to the other ULX pulsars and, indeed, as pointed out in Walton et al. (2018a), very similar to all ULXs in general, regardless of the known or unknown accretor (see also Koliopanos et al. 2017; Pintore et al. 2017). As presented by Walton et al. (2018a), the temperature of the cool component for ULXs ranges from 0.2–0.5 keV. We find for ULX7 that it is 0.33 keV. The hotter component has a temperature range of 1.2–3 keV, and is 2.5 keV for ULX7, making ULX7 indistinguishable from other ULXs in terms of its disk temperatures. Walton et al. (2018a) also noted that the temperature ratio for ULX pulsars in their sample was ∼3, while the other ULXs had a temperature ratio of ∼8. For the ULX7, the ratio is 7.5, therefore more in line with the other ULXs, rather than the other ULX pulsars. Walton et al. (2018a) also found that the flux ratio of the pulsed component, modeled by cutoffpl to the total flux in the 0.3–40 keV band, was higher for ULX pulsars than for the other ULXs. For M51 ULX7, however, this ratio is 0.15, which is relatively low in comparison to NGC 5907 ULX (0.82 in the high state), and NGC 7793 P13 (0.59).
7.3.The Spectral Evolution
The spectral evolution of ULX7 has been studied before. Both Yoshida et al. (2010) and Earnshaw et al. (2016) fitted the available Chandra and XMM-Newton data at the time with a power-law model. Earnshaw et al. (2016) found that the spectral slope did not change significantly, despite the large changes in LX. This included the very faint states observed with Chandra. Rodríguez Castillo et al. (2020) also investigated the evolution of the spectral parameters in their high-quality XMM-Newton data from 2018, using their disk+blackbody model, also finding limited evidence for spectral variations, as did Gúrpide et al. (2021).
In our investigation of the spectral parameters of ULX7 and how they depend on the superorbital phase of the system, we have found that the flux modulations are primarily driven by changes in the flux of the hotter disk-like component, with the flux of the cooler disk component changing less dramatically. We also see a potential dependence of both the disk temperature and normalization of the cooler disk-like component on superorbital phase, even though the flux of this component does not change significantly with the superorbital phase. However, we do not have enough data to claim a significant relationship.
Since the normalization of the diskbb and diskpbb components are directly related to the inclination of the disk (also the distance to the source and inner disk radius), this potential dependence could be straightforwardly interpreted as disk precession, which has been suggested as the mechanism for the superorbital periodic flux variations seen in ULX pulsars (e.g., Dauser et al. 2017; Fürst et al. 2017; Middleton et al. 2018). Indeed, Fürst et al. (2017) also found evidence that the spectral parameters of the disk-like component of NGC 5907 ULX were dependent on the superorbital phase of the system. In this case it was the hot disk-like component, rather than the cool one we see it in, and it was the radial temperature index, p, rather than the normalization/temperature. Since NGC 5907 ULX is relatively absorbed, the cool disk component is not visible. We note, however, that the models used are purely phenomenological, and any physical interpretations should have this added caveat.
The lack of any strong evolution of NH with superorbital period rules out a warped accretion disk that periodically obscures the X-ray source. While small variations in NH can be seen, much larger changes are needed to produce the >1 magnitude flux modulation.
7.4.The Evolution of the Superorbital Period
Having discovered that the X-ray flux from M51 ULX7 is modulated on a period of 38.2 ± 0.5 days in Brightman et al. (2020b) from 2018–2019, we find here that period is variable, and has shifted to 44.2 ± 0.9 days in 2020–2021. Vasilopoulos et al. (2020a) also explored the Swift/XRT data on ULX7, finding the same 38 day period, and also found tentative evidence for a 49 day period from Swift/XRT data taken in 2011, albeit covering only 1–2 cycles and therefore quite uncertain.
An example of a superorbital flux modulation in a source where the period is variable is SMC X-1, which has a superorbital period of ∼60 days (Gruber & Rothschild 1984), and exhibits recurrent excursions to shorter superorbital periods, which may be (quasi-)periodic themselves (Hu et al. 2019). Unlike M51 ULX7, SMC X-1 shows strong spectral variability across its superorbital phase, with a harder spectrum at low fluxes, suggestive of absorption being the main cause. This absorption appears weaker during the superorbital period excursion. From modeling of the pulse profile with superorbital phase from XMM-Newton and NuSTAR data, Brumback et al. (2020) found that the pulse shape and phase of SMC X-1 are consistent with reprocessed emission from a precessing inner disk. The wind-fed HMXB IGR J16493–4348 also exhibits a variable superorbital period; however, in that case, the amplitude is found to be variable, rather than the period. They suggest this is linked to a variable accretion rate. However, they note that the timing and spectral properties of ULXs show significant differences compared to those observed in wind-fed HMXBs such as IGR J16493–4348.
As discussed above, the evolution of the spectral parameters across superorbital phase support the hypothesis that the flux modulations from M51 ULX7 are caused by precession of the accretion disk. Other theories, such as a warped accretion disk periodically obscuring the X-ray source as suggested for SMC X-1, can be ruled out by the lack of NH variations. It has also been suggested that a third orbiting star could be causing the variations (e.g., Middleton et al. 2018; Rodríguez Castillo et al. 2020). This now appears unlikely since if the orbit of the third star were causing the periodic variations, the orbit would have needed to change significantly in a short time.
Vasilopoulos et al. (2020a) discussed the potential for free precession of the neutron star, which has been invoked to explain long-term periodic changes of isolated neutron stars. However, they note that free precession of the neutron star alone cannot account for the variation in superorbital period, which we find here.
8.Conclusions
In our new NuSTAR and XMM-Newton data on M51 ULX7, we have found that the neutron star powering the source has spun up at a rate of 3 ± 0.5 × 10−10 s s−1 since the previous observations by XMM-Newton, which is similar to that seen in other ULX pulsars. The data also provide the first high-quality broadband spectrum, consisting of two disk-like components, the temperatures of which are indistinguishable from other ULXs, and a high-energy tail. We found that the luminosity of the hotter component drives the superorbital flux modulation seen from the source. Finally, we discovered that the superorbital period varies, and has increased from 38.2 ± 0.5 days in 2018–2019 to 44.2 ± 0.9 days in 2020–2021. This change in period rules out some alternative explanations of the superorbital period.
Table 1.Results from the X-ray Spectral Modeling of Joint XMM-Newton and NuSTAR Data on ULX7
| Parameter | Model | ||
|---|---|---|---|
| diskbb+bbodyrad | diskbb+bbodyrad+cutoffpl | diskbb+diskpbb+cutoffpl | |
| NH/1020cm−2 | 1.3 ± 0.7 | 3.3 ± 0.7 |  |
| Tin, cool/keV | 0.46 ± 0.02 | 0.36 ± 0.02 | 0.33 ± 0.03 |
| Ncool/10−4 a |  |  |  |
| Thot/keV | 1.64 ± 0.05 | 1.39 ± 0.09 |  |
| p | ⋯ | ⋯ |  |
| Nhot/10−4 a |  |  |  |
| FX/10−13 erg cm−2 s−1 b | 8.3 ± 0.3 | 8.2 ± 0.3 | 7.7 ± 0.2 |
| LX/1039 erg s−1 c | 7.7 ± 0.4 | 7.6 ± 0.4 | 7.8 ± 0.4 |
| Cpn |  | 0.99 ± 0.07 |  |
| CMOS1 |  |  |  |
| CMOS2 |  |  |  |
| CFPMB | 1.00 ± 0.08 |  | 1.00 ± 0.08 |
| C-statistic/d.o.f. | 4027.81/4286 | 3954.59/4285 | 3951.17/4284 |
Notes. + u indicates a parameter has hit its upper bound in the fit.
a In units of (Rin/D 10)2 cosθ. b Observed in the 0.3–10 keV band. c Intrinsic, corrected for absorption in the 0.2–20 keV band assuming a distance of 8.58 Mpc.
Download table as: ASCIITypeset image
Table 2.Details of the Observations Used in this Work
| Observatory | ObsID | Start Time | Exposure a | Count Rate | Flux b | Superorbital Phase |
|---|---|---|---|---|---|---|
| (UT) | (ks) | (counts s−1) | (erg cm−2 s−1) | |||
| XMM-Newton | 0824450901 | 2018-05-13 21:18:47 | 60.6 | 0.167 ± 0.002 | 5.8 × 10−13 | 0.15 |
| XMM-Newton | 0830191401 | 2018-05-25 20:26:58 | 77.6 | 0.018 ± 0.001 | 2.9 × 10−14 | −0.53 |
| XMM-Newton | 0830191501 | 2018-06-13 01:39:03 | 48.2 | 0.221 ± 0.002 | 7.9 × 10−13 | −0.04 |
| XMM-Newton | 0830191601 | 2018-06-15 01:24:21 | 49.3 | 0.209 ± 0.002 | 7.5 × 10−13 | 0.01 |
| NuSTAR | 60501023002 | 2019-07-10 05:56:09 | 162 | 0.006 ± 0.0002 | 7.6 × 10−13 | 0.26 |
| XMM-Newton | 0852030101 | 2019-07-11 10:47:26 | 58.8 | 0.209 ± 0.002 | 7.6 × 10−13 | 0.26 |
Notes.
a After filtering. b 0.3–10 keV, observed (absorbed).
Download table as: ASCIITypeset image
We thank the referee for their thorough review of our paper which improved it.
This work was also supported under NASA Contract No. NNG08FD60C. NuSTAR is a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS), jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA).
This work was also based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.
We wish to thank the Swift P.I., Brad Cenko, for approving the target of opportunity requests we made to observe M51, as well as the rest of the Swift team for carrying them out. We also acknowledge the use of public data from the Swift data archive. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.
This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC.
G.L.I. acknowledges funding from the Italian MIUR PRIN grant No. 2017LJ39LM. The work of D.S. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. D.J.W. acknowledges support from an STFC Ernest Rutherford Fellowship.
Facilities: Swift(XRT) - Swift Gamma-Ray Burst Mission, NuSTAR - , XMM-Newton. -
Software: NuSTARDAS, XMMSAS (Gabriel et al. 2004), XSPEC (Arnaud 1996), hendrics (Bachetti 2015).
