The variability behavior of NGC 925 ULX-3

We report the results of a 2019-2021 monitoring campaign with Swift and associated target-of-opportunity observations with XMM-Newton and NuSTAR, examining the spectral and timing behavior of the highly variable ultraluminous X-ray source (ULX) NGC 925 ULX-3. We find that the source exhibits a 127-128 day periodicity, with fluxes typically ranging from 1e-13 to 8e-13 ergs/s/cm2. We do not find strong evidence for a change in period over the time that NGC 925 ULX-3 has been observed, although the source may have been in a much lower flux state when first observed with Chandra in 2005. We do not detect pulsations, and we place an upper limit on the pulsed fraction of ~40% in the XMM-Newton band, consistent with some previous pulsation detections at low energies in other ULXs. The source exhibits a typical ULX spectrum that turns over in the NuSTAR band and can be fitted using two thermal components. These components have a high temperature ratio that may indicate the lack of extreme inner disk truncation by a magnetar-level magnetic field. We examine the implications for a number of different models for superorbital periods in ULXs, finding that a neutron star with a magnetic field of ~10^12 G may be plausible for this source. The future detection of pulsations from this source would allow for the further testing and constraining of such models.


INTRODUCTION
Ultraluminous X-ray sources (ULXs) are non-nuclear X-ray sources with luminosities above ∼10 39 erg s −1 , widely agreed at this time to be a population dominated by stellar-mass compact objects accreting at super-Eddington rates (for recent reviews, see Kaaret et al. 2017;Fabrika et al. 2021). Many ULXs are persistently bright but a subset show extreme variability, sometimes of over an order of magnitude in flux, in which they may even enter and leave the super-Eddington luminosity regime. Some of these sources appear to be one-off excursions into the ULX regime, such as particularly bright instances of classical outbursts (e.g. Middleton et al. 2013), or other hard and bright transient ULXs with less certain interpretations (e.g. Earnshaw et al. 2019;Earley et al. 2021). However, many spend significant amounts of their duty cycle in the ULX luminosity regime, or exhibit high-amplitude variability while remaining at ULX luminosities throughout.
One proposed mechanism for such extreme variability in ULXs is the propeller effect (e.g. Tsygankov et al. 2016 for M82 X-2, but see below), in which accretion is halted when the magnetospheric radius of the accreting neutron star exceeds the corotation radius of the accretion disk (Illarionov & Sunyaev 1975;Stella et al. 1986). This would cause the flux to drop dramatically but not necessarily regularly.
Another feature of high-amplitude variability found in some ULXs is (quasi-)periodic variability with timescales of tens to hundreds of days (e.g. Strohmayer 2009; Walton et al. 2016;Vasilopoulos et al. 2020). In some cases, this is distinct from a known orbital period and attributed to superorbital periodicity. Detection of this kind of long-term periodicity requires monitoring over long timescales, and sometimes can be found to be the cause of variability that at first looked like the propeller effect (e.g. Brightman et al. 2019 on M82 X-2). Additionally, some sources demonstrate both long-term periodicity and dramatic flux drop-outs, suggesting that both processes can potentially contribute to the longterm variability of a single source (e.g. Walton et al. 2015Walton et al. , 2016Israel et al. 2017;Fürst et al. 2017). In cases where the variation is confirmed to be superperiodic, the cause of the variability is not yet fully understood, but may be related to the precession within the system, and a variety of models have been proposed to explain it (e.g. Mushtukov et al. 2017;Middleton et al. 2018Middleton et al. , 2019Vasilopoulos et al. 2020).
Interestingly, many known ULX pulsars discovered to date show some form of high-amplitude flux variability, and searching for such variability in archival X-ray data may be a way to identify further ULX pulsar candidates (Earnshaw et al. 2018;Song et al. 2020). However, verification of the nature of this variability requires regular monitoring over a long baseline in order to establish the magnitude of the variability and the presence of any periodicity.
We identified NGC 925 ULX-3 as a highly variable ULX in Earnshaw et al. (2020), henceforth E20. ULX-3 is located at 02 h 27 m 20 s .18, +33 • 34 12.84 (J2000), and was serendipitously detected as a ULX in a Chandra observation of the spiral galaxy NGC 925 at 9.56 Mpc. Upon searching through archival Swift data, we found that it had been bright on a previous occasion, but undetected or detected at significantly lower luminosities since then. Since the sparse existing temporal coverage of observations of NGC 925 made it hard to determine the nature of ULX-3's variability, we undertook a monitoring campaign using Swift, supplemented by target-ofopportunity observations with XMM-Newton and NuS-TAR triggered when ULX-3 was bright.
In this paper, we describe our observation campaign and the reduction of the X-ray data in Section 2, and the results of our analysis in Section 3, including the discovery of a 127-128-day periodicity. We note that a similar result was recently presented by Salvaggio et al. (2022) during the preparation of this paper. Since our analysis is slightly different and we include observations not included in Salvaggio et al. (2022), we believe that this study provides a useful independent confirmation of the periodicity, as well as further spectral analysis. Finally, we present our discussion of the results in Section 4 and conclusions in Section 5. Additionally, all archival Swift data for NGC 925 was used (observation IDs: 00045596001-018), as well as additional monitoring of NGC 925 undertaken between 2021 June 26 and 2021 December 4 (observation IDs: 00014387001-024). Results from our previous analysis of archival XMM-Newton and Chandra data were also used, as detailed in E20. We show the observations used in Table 1.

Swift
We downloaded all Swift data (98 observations in total) and created clean event lists using xrtpipeline. We extracted source and background products with the xselect task, using a 30 radius circular source region centered on the Chandra location of the source (E20) and a 70 radius circular background region located outside the galaxy. Auxiliary response files were created using the task xrtmkarf and the relevant redistribution matrix obtained from the CALDB. Fluxes were calculated from the background-subtracted count rate using PIMMS. As in E20, we assumed an absorbed power-law model, although this time we assumed N H = 1.4 × 10 21 cm −2 , based on our fits to the bright XMM-Newton detection (see Section 3.2.2), since just a The good exposure time after any filtering has been applied, in ks. Given as a range for observation ranges, and as a total for co-added observations. b XMM-Newton exposure times are given for the EPIC-MOS1/EPIC-MOS2/EPIC-pn instruments respectively, and NuSTAR times for FPMA/FPMB respectively. c Made up of Swift observations 00045596035-37 co-added with 00089002001 (see Section 3.2.1).
using the Galactic value would be an underestimate. N H values on the order of 10 21 cm −2 are typical for ULXs (Winter et al. 2006). For the power-law slope we assumed Γ = 1.7 (consistent within the uncertainties for both the high-flux Chandra and low-flux XMM-Newton observations in E20, as well as with the bright XMM-Newton detection fitted by itself). Errors on the count rate (and resulting flux) were calculated in the same manner as in Evans et al. (2009) -that is, for instances of <15 counts, the Bayesian approach of Kraft et al. (1991) was used to calculate the 3σ confidence intervals. If the source was detected, the same method was used to calculate 1σ error bars, otherwise the 3σ upper limit was used.

XMM-Newton
We extracted the data from XMM-Newton observation 0862760201 from the EPIC-pn and EPIC-MOS instruments using the XMM-Newton SAS v18.0.0 software. We produced calibrated event lists using the tasks emproc and epproc. Periods of high background flaring were removed by filtering out intervals of time during which the >10 keV count rate exceeded 0.35 cts/s across the EPIC-MOS detectors and the 10-12 keV count rate exceeded 0.5 cts/s across the EPIC-pn detector. There were no major periods of background flaring during this observation, so this filtering had minimal effect. We extracted data products from a 30 radius circular source region, using a 45 radius circular region on the same chip with a similar distance from the readout node for the background. Events with FLAG==0 && PATTERN<4 were selected from EPIC-pn, and PATTERN<12 from EPIC-MOS. The EPIC-pn light curve was extracted with a bin size equal to the pn time resolution, i.e. 73.4 ms. The spectrum was grouped into 20 counts per bin to allow for Gaussian statistics when fitting, and oversampling set to three times the spectral resolution. The redistribution matrices and auxiliary response files were created using the tasks rmfgen and arfgen, respectively.
The 0.3-10 keV flux for each observation was calculated from the best-fitting absorbed power-law model to the XMM-Newton data by itself, using N H = 1.4 × 10 21 cm −2 as the fixed absorption for the low-flux archival observation.

NuSTAR
The two NuSTAR observations were reduced using NuSTARDAS v2.0.0, with CALDB version 20211020. The nupipeline routine was used to produce clean event files, and nuproducts used to extract source and background spectra and response files. source region at the Chandra source location and centered on the source PSF was used. This is smaller than often typical for NuSTAR extraction regions to minimize contamination from the nearby source ULX-2, which is ∼50 from ULX-3. The background spectrum was extracted from a 60 region located nearby on the chip but outside of the galaxy. Spectra were grouped into 20 counts per bin, as for XMM-Newton.

RESULTS
We show the long-term light curve of ULX-3, both over the course of all X-ray observations and zoomed in on our monitoring campaign, in Fig. 1. The light curve clearly shows repeated brightening and dimming, with the Swift monitoring campaigns covering four complete excursions into a bright state, plus the latter half of a brightening at the beginning of the 2019 monitoring. Each bright state reaches a reasonably consistent peak flux around 8 × 10 −13 ergs s −1 cm −2 (L ∼ 9 × 10 39 ergs s −1 ), with the highest flux we detect dur-ing the monitoring being 1.2 ± 0.3 × 10 −12 ergs s −1 cm −2 (L = 1.3 ± 0.3 × 10 40 ergs s −1 ). Outside of these bright excursions, the source is not always reliably detected by Swift, but upper limits or detections tend to be consistent with a flux of ∼1×10 −13 ergs s −1 cm −2 or a little lower. The lowest-flux detection was by XMM-Newton in 2017, at a flux of 3.5 ± 0.5 × 10 −14 ergs s −1 cm −2 (L = 3.8 ± 0.5 × 10 38 ergs s −1 ). Even lower than this was an upper limit to the flux found by Chandra in 2005 at 2.8 × 10 −14 ergs s −1 cm −2 . The full extent of the observed variability amplitude between the highestflux detection and the lowest upper limit is a factor of 40.

Timing Analysis
To test for long-term periodicity in the light curve, we performed two tests: the Lomb-Scargle test (Lomb 1976;Scargle 1982), as implemented by astropy, and often used for finding signals in sparse or irregularly-sampled data; and an epoch-folding test (Leahy et al. 1983), for which we calculate the L-statistic (Davies 1990). For the epoch-folding test, we folded the flux light curve over each period in eight phase bins, to minimize instances of under-full phase bins.
For each method, we tested for periods between 50 and 200 days, in intervals of 0.1 days, for 1500 trial periods in total. A lower limit of 50 days was chosen because while some ULXs show periodicity on a faster timescale than this, we found it was common for ULX-3 to remain bright for ∼50 days at a time, so we set this as a lower bound. An upper bound of 200 days was chosen since there are at least four excursions into a high-luminosity regime in the ∼800 days of regular Swift monitoring. We included both detections and upper limits from all soft X-ray observatories (i.e. Swift, XMM-Newton, and Chandra) in this analysis, treating upper limits as zeroflux detections. (We obtain the same results at very similar significance if we instead use the measured aperture flux for the observations with upper limits. If we use the method in Salvaggio et al. 2022 and assign each upper limit a flux between zero and the value of the upper limit, the significance of our detection is much reduced, but we found that this effect is mainly due to two observations with weakly constraining upper limits higher than the maximum flux measured for this source. Picking a random flux in the given range for these observations is unrealistic and distorts the results. If these two upper limits are discounted, the significance becomes once again very similar to the zero-flux or aperture-flux treatments of the upper limits).
We performed these analyses for data points only in the 2019-2021 period in which regular monitoring was undertaken, since it is both possible for super-orbital periodicity to change (e.g. Trowbridge et al. 2007;Brightman et al. 2022) or disappear (e.g. Grisé et al. 2013) and for superorbital periodicity to remain consistent even through periods when the flux drops dramatically for other reasons (e.g. Fürst et al. 2017). In each case, the significance of a periodicity was estimated by performing Monte-Carlo simulations following Walton et al. (2016), in which we simulated 10,000 light curves based on a red noise power spectrum at 20 times the total duration of observations and 2 ks resolution (a typical duration for a Swift observation). These light curves were then sampled at a similar distribution of times to the observations, with a scatter of ±1 day applied, and the Lomb-Scargle and epoch-folding tests performed on each. We plot the test statistic results by period, as well as the 4/3σ level where 99.994%/99.7% of the simulated test statistics lie beneath this threshold (Fig. 2, left). As a check, we also performed this analysis on the detections by themselves, leaving out upper limits. We find that we recover approximately the same period, albeit at a lower significance due to the lower number of data points (Fig. 2, right).
To find the error on the resulting periods, we next simulated 1,000 light curves in a similar fashion as above, but this time using as a model the best-fitting sinusoid to the data, with the period set to that found by each test. They were sampled in the same manner and both tests used to find the period for each light curve. We used the distribution of the results from these simulations to determine the 90% confidence interval.
We find that the Lomb-Scargle and epoch-folding tests give consistent results, of 127.2±0.4 and 128±1 days respectively. We show phase-folded plots for both of these periods in Fig. 3. These results are consistent with and confirm the detection of a 126 ± 2 day period first reported by Salvaggio et al. (2022). If we extrapolate these periods back to past data, we find that most previous data points are roughly consistent with such a period, except for the first Chandra non-detection in the case of a slightly longer period (Fig. 4).
We also searched for pulsations in both the XMM-Newton and NuSTAR data using HENaccelsearch from the HENDRICS v5 software package (Bachetti 2018),  searching in the 0.01-5 Hz frequency range and using an accelerated search since ULX pulsars tend to have such dramatic spin-up that a non-accelerated periodicity search would fail to detect pulsations. We did not find any significant detections or potential candidates. Since no tentative signals were apparent from an accelerated pulsation search, we did not perform any more complex analysis involving orbital period modulations, although we note that if the 127-128-day periodicity is an orbital period, we would not expect it to have a dramatic effect on the pulse period modulation over the course of our observation in addition to spin-up. See Section 4.1 for further discussion of the long-term periodicity.
We used the stingray software package to simulate observed light curves of a pulsation at various pulse fractions, following the method described in Fürst et al. (2021), and found that we can place an upper limit on the pulsed fraction of ∼40% in the XMM-Newton band for pulse periods from 0.01 Hz to ∼2 Hz, with the upper limit quickly increasing to 100% for higher pulse frequencies approaching the EPIC-pn timing resolution. There are insufficient data to place a meaningful limit on the pulsed fraction in the NuSTAR band.

Spectral Analysis
We used XSPEC v12.10 (Arnaud 1996) to perform all spectral fitting, and all quoted models are given in XSPEC syntax. In all cases, spectra are grouped into at least 20 counts per bin to allow for χ 2 statistics to be used in fitting. We give uncertainties at the 90% confidence level, and we use the abundance tables of Wilms et al. (2000) throughout. We use only a single tbabs model to account for absorption due to the interstellar medium, though we note that the Galactic contribution to this is N H = 7.26×10 20 cm −2 (Willingale et al. 2013). We considered Swift and XMM-Newton data in the energy range 0.3-10 keV, and NuSTAR data in the range 3-20 keV, above which the background dominates the spectrum for this source. We show all spectral fitting results in Table 2.

NuSTAR Epoch 1
The first NuSTAR DDT observation was taken when the source was not visible to XMM-Newton, so the lowenergy coverage of the spectrum is provided by Swift. In order to increase the signal at low energies, we used the FTOOLS routine addascaspec to co-add the observation taken simultaneously with the NuSTAR observation (observation ID: 00089002001) with the three preceding Swift observations which were all consistent in flux (observation IDs: 00045596035-37). Since there are insufficient low-energy data to place good constraints on N H , we froze it to N H = 1.4 × 10 21 cm −2 , based on our fits to the bright XMM-Newton detection (see Section 3.2.2).
We fitted the spectrum with a number of absorbed single-component models (see Fig. 5 for the best-fitting spectrum and residuals). While a power-law model formally provides an acceptable fit, we find that a cutoff power-law model with a characteristic cut-off energy of ∼6 keV offers a statistically significant improvement (∆χ 2 = 19 for 1 degree of freedom). We can statistically rule out a standard disk blackbody model, though a broadened disk model provides an equally acceptable fit to the cut-off power-law. Using the best-fitting model, we find an absorbed source flux of 5.0 ± 0.4 × 10 −13 ergs s −1 cm −2 (L = 5.5 ± 0.4 × 10 39 ergs s −1 ) in the 0.3-10 keV band. For the purposes of comparison, we also fit the spectrum with two disk blackbody models, although we are unable to place strong constraints on the model parameters.   a The model names are abbreviated as follows: pl = powerlaw, cpl = cutoffpl, dbb = diskbb, dpbb = diskpbb.

NuSTAR
b Value frozen due to the low-energy data being of insufficient quality to constrain NH.
We began by fitting the XMM-Newton spectrum by itself with an absorbed power-law in order to characterize the low-energy emission and quantify the interstellar absorption for use when analyzing the Swift data. We found an absorption of N H = 1.4 × 10 21 cm 2 . While it is possible for N H to change between observations, it is nonetheless likely to be a better approximation to the true N H than the Galactic value by itself. We found that a simple absorbed power-law model offers a statistically acceptable fit to the low-energy data, although there is an 'm'-shaped curvature to the residuals (often seen in ULXs when the low-energy data is fitted with a power-law) that suggests the contribution of multiple components.
We next fitted the XMM-Newton and NuSTAR data simultaneously (see Fig. 6 for the spectrum and residuals). We included a multiplicative constant in this model fit to account for calibration differences between XMM-Newton and NuSTAR, freezing the value to 1 for XMM-Newton and letting it vary for NuSTAR. This constant is anomalously low, with XMM-Newton and NuSTAR expected to be consistent to within ∼10%. The relative normalization may have been affected by the source being very close to one of the NuSTAR chip gaps.
We found that a cut-off power-law offers a significant improvement in fit over a simple power-law, with a slightly higher cut-off energy although still within the uncertainties of that found for the first NuSTAR epoch. A broadened disk model also offers a statistically acceptable fit. However, in all of these single-component cases, the residuals show a characteristic 'm'-shaped curvature indicating that there are likely two components to the emission instead of just one.
We used two-thermal-component models to fit the data, first using two multicolor disk blackbody models, then replacing the higher energy model with a pfree broadened disk model, as is often required for ULX spectra (e.g. Walton et al. 2018). We found that both models fit the data well, although using a broadened disk for the hot component over a standard multicolor  Figure 6. The unfolded XMM-Newton MOS1, MOS2, and pn spectra (black, red, and green respectively) and NuSTAR FPMA and FPMB (blue and cyan respectively) spectra fitted with the const*tbabs*(diskbb+diskbb) model, and the residuals for all fitted models. disk does not show strong evidence for significant broadening and does not provide a statistically significant improvement to the fit. Using the diskbb+diskbb model, we calculated an absorbed source flux of 5.5 ± 0.1 × 10 −13 ergs s −1 cm −2 (L = 6.0 ± 0.1 × 10 39 ergs s −1 ) in the 0.3-10 keV band.

Low Luminosity Swift data
In order to investigate the spectrum of the source at lower fluxes, we used addascaspec to co-add all Swift observations with detections below a flux of 2 × 10 −13 ergs s −1 cm −2 , or with upper limits, discounting low-flux observations or non-detections that are evidently still within an overall bright state (such as the single anomalously low-flux Swift observation in the middle of the most recent bright state). In total, we added 60 observations, covering approximately half of the phase cycle (see Fig. 3). As with the previous Swift analysis, we froze N H to the value measured for the XMM-Newton observation, then we fitted the spectrum with a power-law model. We found that the spectrum could be fitted with a photon index of Γ = 1.8 ± 0.3, with an absorbed flux of 6 ± 2 × 10 −14 ergs s −1 cm −2 (L = 7 ± 2 × 10 38 ergs s −1 ).

A long-term periodicity
We confirmed the high levels of variability of NGC 925 ULX-3 first reported in E20, and discovered a 127-128day periodicity over which the source enters and leaves the ULX luminosity regime. This period is consistent within error with that also reported by Salvaggio et al. (2022). A small number of ULXs, including several of those that have been identified as neutron star accretors, exhibit (likely) superorbital periodicity on the order of tens to hundreds of days (e.g. Strohmayer 2009;Lin et al. 2015;Walton et al. 2016;Brightman et al. 2019), although we note that NGC 7793 P13 has been suggested to have a ∼1500-day superorbital period (Motch et al. 2014;Fürst et al. 2018). We also note that Galactic source SS 433, likely also a super-Eddington accreting source viewed at high inclination, exhibits a 164-day periodicity (Abell & Margon 1979). In this context, the NGC 925 ULX-3 periodicity is fairly typical of the ULX population exhibiting such periods. The flux variation of around an order of magnitude is also within the range of such ULX periodicity discovered so far.
In a couple of instances, the source exhibits a dip in luminosity in the middle of an otherwise bright state, the most obvious of these being about halfway through the most recent bright period, during which there is a Swift detection consistent in flux with measurements made during low-flux intervals. The phase-folded light curves indicate that this dip may recur at a similar phase of the cycle, and likely does not last longer than a few days, although higher-cadence observations would be required to confirm this. This dipping behavior has also been seen in other ULXs with long periods (e.g. Pasham & Strohmayer 2013;Walton et al. 2016), and may be due to periodic/superperiodic obscuration of the source. However, the spectrum of NGC 925 ULX-3 indicates that it is unlikely that we are viewing the source at a high inclination (see Section 4.2).
This period is more-or-less consistent with archival data points, with the main outlier being the Chandra non-detection in 2005. This may be due to a change or disappearance of the super-orbital period in the intervening time, as has been observed elsewhere (e.g. Grisé et al. 2013;Lin et al. 2015;Weng & Feng 2018;Brightman et al. 2022). However, we note that the Chandra upper limit is 2.8 × 10 −14 ergs s −1 cm −2 , significantly lower than the low-flux detections or our composite lowflux spectrum. This may instead indicate an interval of lower flux by some other mechanism, over which a superorbital period may still persist (e.g. Fürst et al. 2021). Therefore, while Salvaggio et al. (2022) claim that the period or phase may have changed in archival observations compared with the periodicity detected in our monitoring campaign, we believe that there is insufficient archival evidence to make such a claim.
There have been various models proposed to explain superorbital periods in ULXs. While an in-depth theoretical exploration is beyond the scope of this paper, we investigated the implications of a 127-128-day period using some of these models and the assumptions within the corresponding papers. Mushtukov et al. (2017) propose that such variability may be caused by superorbital precession of the magnetic dipole of the accreting neutron star, in the context of a system in which the hot thermal emission originates from an accretion curtain that envelops the entire magnetosphere of the ULX, and the cooler thermal emission comes from a supercritical accretion disk. We find that, if we assume a 1-10 s pulsation period typical of ULX pulsars so far discovered, a very high magnetic field strength of >10 15 G is required to produce a superorbital period in the region of 127-128 days, and the corresponding expected temperatures of the thermal components are far lower than we observe (∼0.01 keV). If we assume a pulsation period on the order of minutes rather than seconds, the model requires lower values of the magnetic field strength of 10 11 -10 12 G, which gives more reasonable thermal component temperatures in the region of ∼1 keV, although the two model components are always much closer in temperature to each other than what we observe.
Another proposed model for superorbital variability in ULXs is that of Lense-Thirring precession of the outflowing wind (Middleton et al. 2018(Middleton et al. , 2019, which can be used to construct a timing-accretion plane that may indicate whether the accreting object is a candidate black hole. Here we use the model of a supercritical inner disk and a cooler outer disk, with the observed temperature of the cooler component corresponding to the temperature at the spherization radius T sph . In this case, a precession period of 127 days and T sph = 0.35 keV place the source comfortably within the region of the timingaccretion plane that contains both black hole and neu-tron star accretors (see Fig. 4 of Middleton et al. 2019). We find that, using this model, a precession period of 127-128 days can be generated with a magnetic field strength in the range 10 11 -10 12 G, with the energy fraction used to launch the wind wind = 0.3-0.35. Vasilopoulos et al. (2020) suggest that superorbital periodicity may instead be due to free precession of the neutron star itself as described in Jones & Andersson (2001), in which the distortion of the neutron star can be derived from its spin and precession period and related to the surface magnetic field, accounting for superconductivity in the neutron star interior (Lander 2014), with the precession of the accretion disk synchronized to the neutron star via some coupling with the magnetospheric field lines. If we once again assume a spin period of 1-10 s, a magnetic field of 10 12 -10 13 G is required for this model to reproduce long-term periodicity on the timescale we measure. Weaker magnetic fields would correspond to shorter spin periods.
A detection of pulsations from this source, allowing the determination of the neutron star spin and any spin-up, would help to further test and constrain these various models. We note, however, that these models assume that the magnetic field is dominated by a dipole component -it has been suggested that a significant multipolar component to the magnetic field may be present in at least some ULXs (Israel et al. 2017;Tsygankov et al. 2017), which may introduce further complexity to these scenarios. We also note that the precession of a radiation-driven warped disk, as seen in some other X-ray binaries (Ogilvie & Dubus 2001), has also been proposed as a mechanism for producing a superorbital period for some ULXs (e.g. Kong et al. 2016).
Since an orbital period has not been identified for this source, it is possible that this periodicity is instead orbital rather than superorbital. Several of the ULXs with confirmed orbital periods have periods on the order of days (e.g. Bachetti et al. 2014;Israel et al. 2017), though NGC 7793 P13 has an orbital period of ∼65 days (Fürst et al. 2021), so a period on this timescale is not out of the question. Orbital periods on the order of ∼100 days have been observed in some Be X-ray binaries (BeXRBs; Walter et al. 2015), with luminosities in the ULX regime for some of the brightest outbursts of BeXRBs, and very large variation in flux being common. However, even during the lowest-flux parts of its phase cycle, NGC 925 ULX-3 has luminosity >10 38 ergs s −1 , far more luminous even than most BeXRB outbursts and certainly more luminous than the intervals between outbursts, for which luminosities of <10 36 ergs s −1 are expected. Additionally, BeXRBs tend to have shorter duty cycles due to their eccentric orbits than what we observe for NGC 925 ULX-3, which is close to 50%. Therefore, it seems more likely that the periodicity we observe is superorbital in nature.

Spectral behavior
The spectrum at higher fluxes is typical of a ULX in the super-Eddington ultraluminous state (Gladstone et al. 2009), with a turnover in the NuSTAR band and a complex shape that can be fitted using two disk blackbody components. The fact that we see a high-energy emission component indicates that we are unlikely to be viewing the source at a high inclination -super-Eddington systems at high inclinations tend to show very soft spectra dominated by the cool thermal component due to the hotter central region being obscured (e.g. Urquhart & Soria 2016). The broadband spectra of ULXs in the ultraluminous state will often show broadening in the hotter component, associated with an advection-dominated supercritical accretion disk, as well as a steep power-law excess at high energies (e.g. Bachetti et al. 2013;Mukherjee et al. 2015;Rana et al. 2015) which may be due to emission from the accretion column of a neutron star (e.g. Walton et al. 2018). We do not find any strong evidence of broadening or the power-law tail for NGC 925 ULX-3, although this is likely due to the limited data quality of our NuSTAR observations, for which we do not have many data points above 10 keV.
The spectra between the two epochs are quite consistent in shape and flux, so it is reasonable to assume that the source is in the same state during the first NuSTAR epoch as it is in the second -fitting the first epoch with a two-thermal-component model suggests that the cooler of the two thermal components may be hotter than that of the second epoch, though the quality of the Swift data is insufficient to draw any strong conclusions from this. Also, since the first NuSTAR epoch is well-fitted with a diskpbb model alone, with similar parameters to the hot component of the second epoch, it may be the case that a cool component simply increased in normalization between the two observations.
A temperature for the cooler component of kT ≈ 0.35 keV is fairly typical for the ULX population with good broadband data (0.8 keV is unusually hot for this component, but our measurement for the first NuSTAR epoch is consistent within measurement uncertainties with more typical cooler values). For the hotter component, kT = 2.4-2.8 keV is also reasonable, if on the high end, compared with other known ULXs. If this component originates from a supercritical accretion disk between the spherization radius r sph and an inner radius truncated at the magnetospheric radius r m of an accret-ing neutron star, in the context of the model proposed in Walton et al. (2018), hotter temperatures of this component would correspond to a smaller r m (or even an accretion disk that is not truncated at all, should the accretor instead be a black hole). If this were the case, we would expect to see a broadened hot component (in highly-truncated cases with r m close to r sph , even a supercritical disk may appear as a narrow spectral component as the emission would originate from a limited band of radii). Higher-quality NuSTAR data would be required to explore the shape of the hot component in greater detail and confirm the presence or absence of broadening in the hot component.
In this model, the cooler component would either originate from an outer thin disk whose inner edge is at r sph , or from the radiatively-driven outflow within r sph (given an outer disk's requirement to be sub-Eddington, it is more likely for this component to be the latter for this source, since the luminosity of the cool thermal component alone is ∼2.5×10 39 ergs s −1 ). In either case, the ratio of temperatures between the two thermal components can give some rough indication of the relative sizes of r m and r sph . The temperature ratio between the two thermal components of NGC 925 ULX-3 in the second epoch is 7-8, whereas ULX pulsars tend to have lower ratios of ∼3 . A larger temperature ratio would indicate that r m << r sph , with the inner thick disk not severely truncated by the magnetic field and contributing more emission than the pulsed emission coming from the neutron star accretion column. This may contribute to the dilution of pulsations, making them harder to detect.
Another factor that may contribute to the nondetection of pulsations is geometric beaming by a collimating wind with a large scale height compared to r m , in which the amount of photon scattering within the accretion funnel dilutes the intrinsic pulsed signal. This is a scenario that can be ruled out for observed ULX pulsars with high pulsed fractions (e.g. Mushtukov et al. 2021) but conversely may apply to ULXs in which pulsations aren't detected, particularly those which show evidence of a strong outflowing wind component such as in NGC 925 ULX-3. These potential factors indicate that our inability to detect pulsations from this source does not necessarily rule out a neutron star accretor. (We also note that pulsed fractions in the XMM-Newton band of 10-20% have been observed in other ULXs, consistent with our upper limit of 40%, e.g. Fürst et al. 2016;Rodríguez Castillo et al. 2019).
The low-flux spectrum is consistent in shape with the low-flux XMM-Newton observation analyzed in E20 (which was well-fitted with a power-law model with Γ = 1.8 +0.2 −0.1 ) as well as with the high-flux XMM-Newton spectrum analyzed in this work. Therefore we find no particular evidence of spectral change across the profile of the X-ray period, such as that observed in NGC 5907 ULX-1 (Fürst et al. 2017), although deeper observations of the low-flux regime for this source will be required to place stronger constraints on the shape of the low-flux spectrum and search for evidence of spectral change with flux. Since there is little spectral change despite the source moving into and out of the super-Eddington luminosity regime, this indicates that there is unlikely to be an intrinsic change in the accretion state itself to produce this variation. The luminosity of 10 39 ergs s −1 used to define ULXs comes from the assumption of a ∼10 M black hole -all luminosities that we measure for this source are above the Eddington luminosity for a 1.4 M neutron star, so there is no requirement for this source to be changing between super-Eddington and sub-Eddington accretion states.

CONCLUSION
NGC 925 ULX-3 shows the two-component spectrum typical of a super-Eddington accreting system, and demonstrates long-term, likely superorbital variation that is common in ULXs containing neutron star accretors, with a period of 127-128 days. Various different models of superorbital periodicity in ULXs could potentially explain its behavior, with neutron star magnetic fields up to ∼ 10 12 G providing reasonable solutions. While the source luminosity varies by approximately an order of magnitude in flux, there is no indication of a change in accretion state. Therefore we believe this source to be a strong candidate for being a neutron star ULX and a valuable target for follow-up investigations searching for pulsations.
We thank our anonymous referee for useful suggestions to improve this paper. We wish to thank Dan Stern for useful comments on this work. This work was supported by NASA grants 80NSSC20K1496 and 80NSSC21K0873, as well as by NASA contract NNG08FD60C. TPR acknowledges funding from STFC consolidated grant ST/000244/1. We wish to thank the Swift PI, Brad Cenko, for approving the target of opportunity requests we made to monitor NGC 925. We also wish to thank the NuSTAR PI, Fiona Harrison, for approving the DDT request to observe NGC 925 ULX-3 in December 2019. The scientific results reported in this article are based on observations made by XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA; observations by the NuSTAR mission, a project led by the California Institute of Technology, managed by JPL, and funded by NASA; and observations by the Swift mission, as well as public data from the Swift data archive provided by the UK Swift Science Data Centre at the University of Le-icester. We thank the observing teams for these missions for carrying out these observations.