A new candidate pulsating ULX in NGC 7793

We report here the discovery of NGC 7793 ULX-4, a new transient ultraluminous X-ray source (ULX) in NGC 7793, a spiral galaxy already well known for harbouring several ULXs. This new source underwent an outburst in 2012, when it was detected by \textit{XMM-Newton} and the \textit{Swift} X-ray telescope. The outburst reached a peak luminosity of 3.4$\times 10^{39}$ erg\ s$^{-1}$ and lasted for about 8 months, after which the source went below a luminosity of $10^{37}$ erg\ s$^{-1}$; previous \textit{Chandra} observations constrain the low-state luminosity below $\sim$ 2$\times 10^{36}$ erg\ s$^{-1}$, implying a variability of at least a factor 1000. We propose four possible optical counterparts, found in archival HST observations of the galaxy. A pulsation in the \textit{XMM-Newton} signal was found at 2.52 Hz, with a significance of $\sim3.4\,\sigma$, and an associated spin-up of $\dot{f} = 3.5\times10^{-8}$ Hz.s$^{-1}$. NGC 7793 is therefore the first galaxy to host more than one pulsating ULX.


INTRODUCTION
Ultraluminous X-ray sources (ULXs) are extragalactic non-nuclear sources whose luminosity exceed the Eddington limit for a stellar mass compact object (∼ 3 × 10 39 erg.s −1 ), while their non-nuclear position excludes AGNs. They were first discovered using the Einstein Laboratory (Long et al. 1981), the term being coined later by Makishima et al. (2000). To explain such high luminosities without exceeding the Eddington limit, these sources were at first thought to reveal the presence of particularly massive black holes (e.g. Colbert & Mushotzky 1999). However, two discoveries in the following years challenged this interpretation of ultraluminous sources. The first NuSTAR observations of ULXs in the hard X-rays allowed to detect a high-energy cutoff, which was thought to indicate super-Eddington accretion (Bachetti et al. 2013;Walton et al. 2013). Then, the discovery of coherent pulsations in a ULX in M82 (Bachetti et al. 2014) lead to the conclusion that at least some ULXs could host accreting neutron stars, which are much less massive. This implies that the accretion rate in these sources goes well above the Eddington limit, up to a factor 500 for the pulsating ULX in NGC 5907, for instance (Israel et al. 2017a). Thus, ULXs are precious tools for understanding accretion mechanisms beyond simple disk models (Shakura & Sunyaev 1973). For a thorough review of the link between ULXs and Super Eddington accretion, see for example Kaaret et al. (2017).
To date, only 7 pulsating ULXs have been found: M82 X-2 (Bachetti et al. 2014), NGC7793 P13 (Fürst et al. 2016), NGC5907 ULX-1 (Israel et al. 2017a), NGC300 ULX-1 (Carpano et al. 2018), NGC1313 X-2 (Sathyaprakash et al. 2019), M51 ULX7 (Ro-★ E-mail: erwan.quintin@irap.omp.eu dríguez Castillo et al. 2020). Swift J0243.6+6124 has been proposed as the first pulsating ULX in our own galaxy (Wilson-Hodge et al. 2018). These sources have several common properties: mainly, their extreme luminosities (even among ULXs), and their long-term temporal variability. Indeed, they exhibit two levels of emission, with drops of up to three orders of magnitude in flux between highactivity and low-activity episodes (see for instance Walton et al. 2015). This common feature lead to the use of the long term variability of ULXs as a proxy to find more pulsating ULXs (Earnshaw et al. 2018;Song et al. 2020). Identifying the ULXs with neutron star accretors would allow us to determine the fraction of ULXs that are powered by accretion onto neutron stars.
NGC 7793 P13 is one of those pulsating ULXs. It was discovered during the first X-ray survey of NGC 7793 by Read & Pietsch (1999), using the ROSAT telescope. This spiral galaxy lies in the Sculptor group, at about 3.9 Mpc (Karachentsev et al. 2003). The initial survey identified 27 interesting sources in the vicinity of NGC 7793, P13 being the most luminous among them at the time of the observation in 1992, at about 9 × 10 38 erg s −1 . Among those 27 ROSAT sources, another source was later identified as a highly variable ULX, NGC 7793 P9 (Hu et al. 2018a). Both P9 and P13 lie in the outskirts of the galaxy. P13 has been observed a total of 14 times between 2012 and 2019 with XMM-Newton, for a total of over 600ks of exposure. This ULX was shown to harbour a 0.42s pulsar (Fürst et al. 2016;Israel et al. 2017b), and the companion star has been identified as a B9Ia (Motch et al. 2014).
We have discovered another interesting source in this galaxy, NGC 7793 ULX-4. Thanks to the presence of NGC 7793 P13, there have been numerous observations of this galaxy, and consequently of this source. The low number of detections however indicate high The green circle indicates the newly discovered ULX in its highest detected state. This 40" radius circular region was used to extract the photons in the analysis. The blue circles correspond to the known ULXs in NGC 7793, respectively P9 on the top left and P13 on the bottom right, which was the target of the observation and in a low state at this time. The field of view is about 360"×560", and roughly covers the optical extent of the galaxy. variability; taking into account the peak flux and assuming its association with the galaxy implies a highly variable ULX.
In this paper, we present the study of this new source, using a wealth of multi-wavelength data from XMM-Newton, the Neil Gehrels Swift Observatory (hereby Swift), and Chandra in X-rays and from the HST and the VLT Multi Unit Spectroscopic Explorer (hereby MUSE) in the optical. In Section 2, we describe the data selection and analysis. The timing and spectral analysis results of the study are then shown in Section 3. These results are finally summed up and discussed in Section 4.

XMM Newton
While looking for long term variability of XMM-Newton sources in the 4XMM-DR9 catalog (Webb et al. 2020), we came upon the extra-nuclear source 4XMM J235747.9-323457 (RA=359.4496, DEC=-32.5826). It was detected in a single XMM-Newton observation, ObsID 0693760101 (see Figure 1). Since NGC 7793 has been observed by XMM-Newton a total of 14 times between May 14 th , 2012 and November 22 nd , 2019, the HILIGT upper limit server 1 was employed to determine the flux upper limits in the case of non-detections. We chose a power-law spectral model with a Γ = 1.7 slope, a hydrogen column density of = 3 × 10 −20 cm −2 , and working on the 0.2-12 keV band, the server returned 3 flux upper limits for 9 observations (see Table 1). These values were chosen to be similar to the parameters used to determine the fluxes of sources detected in the 4XMM-DR9 catalog. Four observations returned automatically computed detections, which we did not take into account as they were polluted by noise or readout streaks, and consequently are not part of the 4XMM-DR9 catalog. The server also allowed us to retrieve flux upper limits for the XMM-Newton Slew survey, on 3 Slew observations (all details are in Table 1).
We then studied the single XMM-Newton detection of this source, in the observation 0693760101, in May 2012. We used the 1 http://xmmuls.esac.esa.int/upperlimitserver/ XMM Newton Science Analysis System v.18.0.0 2 to extract the data from all three EPIC instruments using epproc and emproc, and we applied the following process for the EPIC-pn data (and for the EPIC-MOS data respectively). We removed the high-flaring periods when the events with pattern 0 and with energies greater than 10 keV yielded a rate greater than 0.4 counts.s −1 (0.35 counts.s −1 for MOS), as is recommended in the XMM-Newton data analysis threads. This left about 29 ks from the initial 50 ks. Events were further filtered by removing the bad pixels and the patterns ≤4 (resp. patterns ≤12). The source spectra were then extracted from a circular region centred on the source with a radius of 40" for all three instruments, in order to reach encircled energy fractions of about 90%. The background spectra were selected from 60" radius circular regions in empty zones on the same CCD as the source. We generated the redistribution matrices and ancillary files using respectively the rmfgen and arfgen tasks. Finally, all spectra were rebinned using grppha, with bins of at least 20 photons, and the fits were performed using a 2 statistic.
For the timing analysis, only the EPIC-pn data was used. No background was subtracted from the data. The data were barycentered using the position of the source and the SAS task BARYCEN.

Swift
A Swift counterpart was found in the 2SXPS catalog , named 2SXPS J235747.8-323458. There are 13 Swift detections, allowing us to study the evolution of its luminosity after the outburst. In 181 additional observations, the source was not detected. The unabsorbed fluxes for each detection were those from the catalog, using a powerlaw with a spectral index of 1.7 and galactic absorption. We used the HILIGT upper limit server to retrieve the flux upper limits for Swift, using the same parameters as for XMM-Newton but in the 0.3 -10.0 keV band. In the same way as for XMM-Newton, we removed the server detections that are not part of the 2SXPS catalog. The observations where the source was detected are described in Table 1.
For the spectral study, we retrieved the source products from the automatic pipeline (Evans et al. 2009). Our goal was to determine whether there was an evolution of the spectral shape between the high state, as observed by XMM-Newton, and the declining state, as observed for some of the Swift detections. Since the Swift observations cover a large part of the outburst, we stacked the detections into two groups: those corresponding to the high state, for which the catalog gave a luminosity above 10 39 erg.s −1 ; and those corresponding to the fainter state, below 10 39 erg.s −1 . For both groups, we regrouped the data so each spectral bin contained at least two photons. The fits were performed using C-statistic instead of 2 , given the low number of photons available.

Chandra
The X-ray observatory Chandra (Weisskopf et al. 2000) has observed NGC 7793 four times, but the source went undetected in all occasions. We therefore estimated 3 upper limits on the unabsorbed flux in the 0.5 -7 keV band using the task srcflux in CIAO version 4.12 with CALDB 4.9.3. As for the Swift data, we assumed a powerlaw of photon index 1.7 subject to neutral absorption along the line of sight (n HGal = 3.38 × 10 20 cm −2 ). The observation details and the resulting upper limits can be found in Table 1.

HST
In order to determine whether our source is indeed located in NGC 7793, we looked for its optical counterpart. We examined three Hubble Space Telescope observations in the B, V and I bands, which were made on December 10 th , 2003 (details are in Table 1). The excellent angular resolution allowed us to look for optical counterparts within the XMM-Newton 1 position error zone. The astropy (Astropy Collaboration et al. 2013 and photutils (Bradley et al. 2020) Python packages were used to compute photometry for the sources within this circle. For each possible counterpart, the source and background fluxes were extracted from circular regions of respective radii 0.08" and 0.16" centred on the counterparts. We assumed a poissonian error on the photon counts.

MUSE
In addition to the different counterpart candidates we retrieved from the HST data, we took advantage of the fact that the MUSE instrument (Bacon et al. 2010) on the Very Large Telescope (VLT) observed the galaxy NGC 7793 during its Adaptive Optics Science Verification run 3 . We used the cube ADP.2017-06-16T14:34:33.584, which had a 3500s exposure time. Using the Python package mpdaf (Bacon et al. 2016), we corrected a 3" offset with the HST data, by assuming a correct HST astrometry and using the estimate_coordinate_offset function, which looks for the peak of the cross-correlation between both images. We then proceeded to extract the optical spectrum (475-935 nm) of the zone associated with the X-ray source, that contains all the HST counterpart candidates. The lower spatial resolution (PSF∼0.75") precludes us from studying individually each of the sources found in the HST image; the spectrum thus corresponds to the combination of the four candidates. Moreover, this source is quite close to the core of NGC 7793, about 40" away, to be compared with 2' for NGC 7793 P13. Considerable extended emission from the galaxy is thus present in this region, making the background subtraction challenging; as a consequence, the spectrum is very noisy. It was not possible to identify the nature of the stellar companion as was done for instance for NGC 7793 P13 (Motch et al. 2014). It was however possible to check for the presence of strongly redshifted or galactic lines, in order to assess the possibility of a foreground or background optical counterpart.

Long term evolution
The long term evolution of the X-ray source is shown in Figure  3 ). The amplitude of variability is about three orders of magnitude, which is consistent with the most variable ULXs in Song et al. (2020).

Short term evolution
In order to analyse the short-term X-ray variability of the source, we focus on the XMM-Newton observation 0693760101, where the source was clearly detected. Figure 4 shows the pn lightcurve corrected using epicclcorr binned to 100s, where small data gaps (t<20s) have been replaced by Poisson distributed rates around the mean count rate computed from the GTIs. As can be seen, the source shows little variability throughout the observation. For the subsequent timing analysis, we selected the events in the time window shown in Figure 4 between the black dashed lines, to avoid artifacts introduced by the large data gaps created by the GTI. This results in a net exposure of ∼ 29 ks.
We next produced two power density spectra (PDS) from the total (i.e. source + background) 29 ks uncorrected lightcurve using the Leahy normalization (Leahy et al. 1983), as we find that the epicclcorr task can alter the underlying statistics. The first one is produced from the entire lightcurve rebinned at the maximal resolution (∼ 74ms), suitable to look for coherent pulsations up to ∼ 7 Hz. For the second one, we bin the lightcurve to 0.5s and average together power spectra of time series of 500s long (see van der Klis (1988) for a review), in order to look for quasi-periodic oscillations similar to those reported by Atapin et al. (2019). We find both PDS consistent with being dominated by Poisson noise. The first PDS is shown in Figure 4 to highlight that no pulsations were identified.

Pulse searches
Due to the high spin up rates commonly found in PULXs, more refined timing techniques are generally required to detect any pulsations. We therefore undertook an accelerated pulsation search using the HENDRICS timing software (Bachetti 2018). We ran an accelerated Fourier-domain search (Ransom et al. 2002), implemented in the tool HENaccelsearch, between 0.01 and 7 Hz. To overcome the loss of sensitivity of the power density spectrum (PDS) close to the borders of the spectral bins, we used "interbinning", a way to interpolate the signal between two bins using the signal in the adjacent bins (van der Klis 1989). This search technique produces many false positives, due to the alteration of the white noise level of the PDS, that need to be confirmed with independent methods. We found 88 candidates, with frequencies between 0.01 and 5 Hz. A few of the ones at the lower boundary of this range turned out to be clear artifacts, and we could reject them by visually inspecting the folded profiles.
We then analyzed the − parameter space around all candidate frequencies with the 2 1 (Rayleigh test) and 2 2 search 4 (Buccheri et al. 1983). To perform this analysis, we used the tool HENzsearch with the fast option, that decreases the computational time of the search by a factor of ∼ 10 and automatically optimizes the search in the space by pre-binning the photons in phase and computing the 2 statistic using these bins instead of the individual photons (Huppenkothen et al. 2019).
As expected due to the use of interbinning, most of the candidates were not confirmed as significant by the 2 search. However, Figure 2. Top Left: MUSE data, averaged. The white circle corresponds to the X-ray position with a 1 error radius of 0.29", and the field of view is roughly 10"×13". Top Right: HST data in the B band. Bottom Left: HST data in the V band Bottom Right: HST data in the I band. All the HST images correspond to a field of view of about 2"×2.5", the white circle still corresponds to the X-ray 1 position error. We can see the four possible counterparts, each with a 0.8" circle around it: Blue circle is Source 1; Red circle is Source 2; Cyan circle is Source 3; Green circle is Source 4. one candidate stood out, with a nominal significance of ∼ 3.4 . The pulse being consistent with a sinusoid, the candidate is detected best in the 2 1 test. We show the details in Figure 5. This candidate has a frequency of 2.5212(1) Hz and a frequency derivative of 3.5(2) · 10 −8 Hz/s. The pulsed fraction is around 12 ± 3% over the XMM-Newton energy band.
We caution that estimating the significance in such a search is tricky. In particular, it is difficult to estimate the correct number of trials after running an accelerated pulsar search where the spectral bins are not completely independent, using interbinning, and then running the 2 2 search on the most promising candidates. Another word of caution: the candidate pulsation is only a few times the frame time. As an additional test, we randomized the event arrival times adding a time delta uniformly distributed between ±36.5 ms (plus/minus half the frame time of EPIC-pn in full frame mode), and doing so the significance of the detection decreased slightly below 3 . Since we are dealing with a small number of counts, it is hard to say whether we artificially diluted the signal of a bona fide candidate.

X-ray spectrum
For the spectral analysis, we used xspec (Arnaud 1996). The spectra obtained from the three XMM-Newton instruments were simultaneously fitted in the 0.2-12 keV range with an absorbed double multi-color black body (tbabs(diskbb+diskbb)), which is commonly used for ULXs (see for instance Koliopanos et al. 2017, and references therein). The spectrum is shown in Figure 6. The results of the fit are shown in Table 2. The temperature values associated with the two black bodies correspond to typical values found for other ULXs, i.e. a cool (T∼0.2 keV) and a hot (T∼1.5 keV) thermal component (see Koliopanos et al. 2017, for more details). To assess the validity of the use of two thermal components, we tried fitting the XMM-Newton data with a single absorbed diskbb; this lead to a much larger value of 2 (369.61 with 326 degrees of freedom, compared to a 325.32 with 324 degrees of freedom for the double diskbb), with a significant soft excess (residuals>5 below 0.6 keV), thus justifying the use of a second low-energy diskbb. Table 1. Log of observations of NGC 7793 ULX-4. For the sake of readability, we did not add the 171 non-detections from Swift (among which 10 were considered as detections by the upper limit server but are not part of 2SXPS, and were therefore ignored). The unabsorbed XMM-Newton detection flux was computed based on the spectral study we performed. The unabsorbed Swift fluxes were taken from 2SXPS, using a fixed 1.7 powerlaw. All XMM-Newton and Swift upper limits were computed using HILIGT, with a 1.7 powerlaw and = 3 × 10 −20 cm −2 . The upper limits correspond to a 3 significance. All Swift observations marked with a dagger were used for the spectral study of the low-luminosity phase, as for them the pipeline yielded luminosity below 10 39 erg.s −1 . The Chandra upper limits were computed using the CIAO tools. The unabsorbed fluxes in different bands were computed using the model component cflux, and freezing the normalisation of any of the black bodies (both yield the same results). The bands used are 0.2 -12 keV to assess the total flux, and 0.3 -2 keV and 2 -10 keV to quantify the hardness of the emission; the hardness ratio being defined as 2−10 − 0.3−2 / 2−10 + 0.3−2 . These latter bands were chosen in order to be compared with Swift data. From the fit, we estimated a 0.2 -12.0 keV unabsorbed flux value of (1.93±0.05) ×10 −12 erg.s −1 .cm −2 ; assuming a distance of 3.9 Mpc (Karachentsev et al. 2003), this implies a peak luminosity of (3.43 +0.16 −0.12 ) × 10 39 erg.s −1 . The Swift high state and low state spectra were then compared to the XMM-Newton spectrum. For both Swift groups, the spectral range used was 0.3 -7.4 keV, since above this value the source flux dropped below the background. The models were then extrapolated to the complete 0.3 -10.0 keV range. In both cases, the hydrogen column density was set to the value found in the XMM-Newton high state observation, to improve the robustness of the fits on this low-quality data; this value is compatible with the galactic value.
The Swift high state spectrum was fitted with the same absorbed double multi-colored black bodies model as the XMM-Newton data, showing comparable spectral parameters to those for the XMM-Newton data. The large error bars, especially for the low energy black body, are due to the smaller Swift collection area and its higher minimum energy. We estimated the hardness of the spectrum in the  same way as for the XMM-Newton data, and the results are shown in Table 3.
Then, the low state spectrum was studied. The idea behind this study was to assess the presence of a spectral change along with the change of luminosity. To this effect, we looked for the minimal changes needed in the high state spectrum in order to account for the low state spectrum. We began by freezing the model parameters from the XMM-Newton data, and only accounting for the luminosity change by adding a constant offset, which was the only free parameter. This lead to a relatively poor fit, with a C-stat of 68 for 66 degrees of freedom. Then, we removed the constant offset and let both black bodies free, but this lead to poorly constrained parameters, both normalisations being only upper limits; freezing the temperature of the cool black body at 0.1 keV prevented this for the hot black body, but not for the cool one. Finally, seeing the poor constraints on the normalisation of the cool black body, we removed it entirely, allowing to have a good fit statistic with no parameters for which the normalisation could reach zero. This final model still allowed us to evaluate the hardness of the emission, obtained by comparing the fluxes of this model in the 0.3 -2 keV and 2 -10 keV bands. The spectrum during the decline phase was in the end consistent with a spectral change, and revealed a spectral softening of the source with decreasing luminosity (see Table 3).

Optical counterpart
Four candidate counterparts were found in the HST data as shown in Figure 2. Their absorbed magnitudes are shown in Table 4. Sources 1 and 2 have high fluxes in the B band (about 3 × 10 −19 erg.s −1 .cm −2 ), while Source 3 is highest in the I band (about 1.5 × 10 −19 erg.s −1 .cm −2 ); Source 4 has a flatter spectrum. It is impossible at this point to determine precisely which star is the  companion to NGC 7793 ULX-4. In order to get a better assessment of the nature of the source, we focused on the X-ray to optical ratio (see Zolotukhin et al. 2016, and references therein). For ULXs, it should reach between 10 2 and 10 3 , about 10 for AGNs and below 0.1 for stars.
Since there are no simultaneous optical and X-ray observation of the source, we first assume a constant optical flux, between the 2003 HST observation and the 2012 ULX outburst; alternatively, we use a Chandra non-detection three months before the optical observation to provide an upper limit on the X-ray to optical ratio of the source in its low state. The typical timescale of the outburst (about seven months) is larger than the interval between the Chandra upper limit and the HST observation, which justifies the use of this second method as as quasi-simultaneous multi-wavelength observation. For the optical flux, we considered the December 2003  Figure 6. Spectrum obtained from XMM-Newton using all the EPIC data: EPIC-pn (black), EPIC-MOS1 (red) and EPIC-MOS2 (green), binned in groups of at least 20 photons, and then simultaneously fitted with an absorbed double multi-colored black body model. For the plot, data was rebinned to 5 . Table 3. Flux values in different bands, obtained using the model component cflux, and the corresponding hardness ratios. The XMM-Newton data corresponds to the peak, in ObsID 0693760101.

Data
Flux 0.3-2 keV Flux 2-10 keV HR (×10 −13 erg.s −1 ) (×10 −13 erg.s −1 ) V band HST observation, computed for each of the optical counterparts. The resulting approximations of the / opt for the different counterparts and both methods are shown in table 5. The approximation of constant optical flux gives very large values of / opt , consistent with ULXs and excluding AGNs or stars. Even if we conservatively assumed a simultaneous optical variability of a factor 10 3 at the time of the X-ray outburst (which is the order of magnitude of the X-ray variability), the ratio indicates a ULX in the galaxy NGC 7793. The / opt ratio for the low state (using the   Table 5. / opt calculated using two different methods: either assuming a constant optical flux, or using the X-ray flux upper limit quasi-simultaneous to the optical observation (three months difference).

XMM-Newton
/ opt method Source 1 Source 2 Source 3 Source 4 opt constant 8.9 × 10 6 8.3 × 10 6 1.3 × 10 6 1.3 × 10 6 upper limit <2.1 × 10 3 <2.0 × 10 3 <2.5 × 10 3 <2.5 × 10 3 Chandra non-detection) gives an upper limit that is consistent with a ULX in NGC 7793, but does not exclude a background AGN or a foreground star. An additional step in confirming the association of these sources with the galaxy is the study of their redshift. MUSE's spatial resolution only allows us to study these sources as a whole, preventing us from determining each source's redshift individually. However, the redshift of both and lines, determined by fitting gaussian profiles to the spectrum, is compatible with that of NGC 7793, which is (7.57 ± 0.07) × 10 −3 (Lauberts & Valentĳn 1989). This further indicates that the source belongs to NGC 7793. The luminosity is therefore compatible with a ULX and the strong variability supports the identification as a transient ULX.

DISCUSSION & CONCLUSION
The first step in studying this source is confirming its association with NGC 7793 by excluding a background or foreground contamination. Several arguments can be made against a background AGN. Firstly, the amplitude of the X-ray variability, about three orders of magnitude, is unlikely for an AGN (which typically have X-ray variability of the order a factor 10, see McHardy 2010). The X-ray softening during the declining phase, shown in Table 3, is the opposite of what is expected for an AGN, which are typically softer during high-luminosity episodes. Additionally, the peak X-ray to optical ratio is very high for an AGN. Finally, the absence of highly redshifted H or H lines in the MUSE data tends to exclude a background galaxy as well. We then consider a foreground source, either a Galactic X-ray binary or a star. The Galactic latitude of about -77 • means that if this source is galactic, it could not be further than about 1 kpc from Earth. Taking the peak flux measured by XMM-Newton, this yields a peak luminosity below (2.26 ± 0.05) × 10 32 erg.s −1 . This peak value is too low for X-ray binaries, which are expected to reach about 10 38 erg.s −1 (Kuulkers et al. 2003). It would be consistent with a stellar flare, but the duration of the peak (about 300 days) is much longer than what is expected for a star (a few hours at most, Pye et al. 2015). The peak X-ray to optical ratio is once again too high for a star. All these arguments tend to exclude any foreground or background object. In addition to this, the spectral shape of the peak emission as observed by XMM-Newton is consistent with that of a ULX, and the long-term variability is comparable to that of other varying ULXs (Earnshaw et al. 2018;Song et al. 2020). This level of variability on such timescales may be due to a central neutron star transitioning from a propeller state (Tsygankov et al. 2016), during which its magnetic field stops the accretion flow and thus limits the luminosity.
The spectral fits allow us to put some constraints on the source; more precisely, the normalisation of the cool diskbb component gives us information about the inner radius of the disk (in km), through the relation = ( disk / 10kpc ) 2 cos( ), being the inclination (Mitsuda et al. 1984;Makishima et al. 1986). Using this relation for the XMM-Newton data, we get that in the high state disk = 550 ± 320km/cos( ) 1/2 .
Due to the unknown inclination, this is a lower limit on the inner radius of the disk. Such a value is comparable with that of other ULXs (Koliopanos et al. 2017). Assuming that the inner radius of the disk coincides with the magnetospheric radius , we can get an estimate of the central neutron star's magnetic field, using the Equation (1)  = 7 × 10 7 Λ 1/7 10/7 6 4/7 12 −2/7 39 cm, where Λ is a constant depending on the accretion geometry, often taken to be Λ = 0.5 for disk accretion; is the mass of the neutron star in units of solar masses, taken to be = 1.4; 6 = /10 6 cm is the radius of the neutron star, taken to be 6 = 1; 12 = /10 12 G is the surface magnetic field strength, and 39 = /10 39 erg s −1 is the accretion luminosity, taken to be the peak X-ray luminosity of 39 ∼ 3. Assuming cos( )∼1, this leads to an estimate of the surface magnetic field of the neutron star during the peak at ∼ 3.5 +4.1 −2.5 × 10 12 G. This value of the magnetic field is once again comparable to that of other ULXs, among them P13, with P13 ∼ 4.2 × 10 12 G (Koliopanos et al. 2017). For the low state, if we again assume that the disk corresponds to the cool black body we can only use the second model from Table 2, which is the only one with a cool component. Its normalisation, although poorly constrained, is consistent with an increase between the high state and low state; this would mean an increase in the inner disc radius, as the lower luminosity implies that the accretion rate drops. The ram pressure then diminishes, so the disc can not penetrate as close to the neutron star.
Another point that can be made about this source's declining state spectrum is the possible presence of emission lines. Indeed, taking a look at Figure 7 reveals that some data bins collectively diverge (∼ 2 ) from the continuum emission, at ∼ 0.7 keV and ∼ 1.6 keV. Taking into consideration the quality of the data in hand, it is impossible to conclude whether these data spikes are indeed emission lines. However, such emission lines around these energies are expected to be detected in ULXs, generally massively blueshifted (Pinto et al. 2017), thus revealing the presence of a powerful wind at relatively large fractions of .
Once the source was confidently identified as a varying ULX, the rest of our study of NGC 7793 ULX-4 focused on the search for pulsation. As stated above, we find a candidate pulsation at 2.5212(1) Hz with an associated frequency derivative of 3.5(2) · 10 −8 Hz/s. Under the hypothesis that this candidate is a real pulsation, the frequency derivative we found would be indicative of orbital motion. Such an would indeed be unrealistically high for the accretion torque on a neutron star, given the flux we measured in Section 3.3.1. According to standard accretion theory (e.g. Ghosh & Lamb 1979;Wang 1995), given a corotation radius co and a magnetospheric radius co one expects the frequency derivative to be roughly related to luminosity (see, e.g. Bachetti et al. 2020) where is 6 for thin disk accretion à la Shakura & Sunyaev (1973) and 3 for super-Eddington accretion à la King (2009). The PULX with the highest luminosity observed so far, NGC 5907 ULX-1 (10 41 erg s − 1, Israel et al. 2017a), has ≈ 10 −9 Hz/s. The peak luminosity of NGC 7793 ULX-4 being much lower, such a high value for the frequency derivative most likely comes from orbital motion. We do not have enough counts to test the presence of a second derivative, which would allow better constraints on the orbital parameters. The measured is in the range observed for the ULXs with a known orbital period in the range of ∼few days (Bachetti et al. 2014;Rodríguez Castillo et al. 2020).
Like M 82 X-2 (Bachetti et al. 2014), the orbital period of a few days implies an intermediate mass stellar companion for mass transfer to occur. Intermediate mass X-ray binaries have frequently been suggested to be at the origin of ULXs (see e.g. Misra et al. 2020, and references therein). Indeed, for a neutron star ULX with a luminosity about 10 times Eddington, the most likely system is an intermediate mass X-ray binary with a period of the order of days (Misra et al. 2020).
The pulse period of NGC 7793 ULX-4 is also surprisingly close to that of P13 (2.4 Hz, Israel et al. 2017b), but the values are extremely different (2.2 × 10 −10 Hz.s −1 for P13 compared to 3.5 × 10 −8 Hz.s −1 for ULX-4), which excludes a possible contamination of the observation by P13. Finally, the fact that P13 was in a lowluminosity state at the time of the peak luminosity for ULX-4 (see Figure 1), and that the two ULXs are separated by 155", further excludes the possibility of a contamination by P13.
Assuming that the pulsation detected for ULX-4 is confirmed, NGC 7793 is the first known galaxy hosting two pulsating ULXs. Another ULX known in this galaxy, P9 (Hu et al. 2018b) has been proposed to contain a black hole based on spectral analysis. However, the strong variability in this source is reminiscent of the large levels of variability observed in PULXs (Song et al. 2020), frequently ascribed to the transition to the propeller regime. A detailed search for pulsations will be the subject of a future paper. Another object exists in NGC 7793 (P8 in Read & Pietsch 1999) that has been proposed to be a third ULX, S26 (Soria et al. 2010). Whilst it is not currently seen as a ULX using X-ray observations, the power output measured using radio observations of the surrounding environment indicates that this is also a ULX. The radio and X-ray observations support a black hole hypothesis for the accretor. No other sources have been formally identified as ULXs in this galaxy. If there are other ULXs, they are likely to be highly variable, as we have not seen them in decades of observations, again indicative of ULXs containing a neutron star. Such a large proportion of neutron star ULXs stands out from other galaxies which have both persistent and transient ULXs. Why NGC 7793 should be different is unclear. Following population synthesis studies, Wiktorowicz et al. (2019) propose that more neutron star ULXs exist either in 1) galaxies with constant star formation rates, approximately solar metallicities and ages of at least 1 billion years, or 2) in galaxies a long time after a burst of star formation in regions with lower than solar metallicity. The central regions of NGC 7793 (in which lies ULX-4) fulfill the criteria for the first case (Vlajić et al. 2011) and the outer regions of NGC 7793 (in which lies P13) fulfill the criteria for case two (Vlajić et al. 2011). This may explain the propensity for neutron star ULXs in NGC 7793.
Other observations of NGC 7793 ULX-4 during an outburst are needed to confirm the pulse period and determine the first and second period derivatives. Further spectral studies over the outburst would help to confirm the physical nature of the processes behind the evolution of luminosity, and would also allow us to rigorously study the presence emission lines and their possible blueshift. Finally, optical observations during the outburst would be useful to identify which of the possible counterparts is the true counterpart. This source has been quiescent since its first and only detected outburst, in May 2012; the constant monitoring of the nearby ULX NGC 7793 P13 will however undoubtedly allow us to detect a potential new high activity period in ULX-4.