© The Authors 2023

1. Introduction

Ultraluminous X-ray sources (ULXs) – defined as extragalactic, off-nuclear, point-like X-ray sources with luminosities in excess of the Eddington limit for a stellar-mass black hole (BH; L_X ≳ 10³⁹ erg s⁻¹, under the assumption of isotropic emission) – have been explained as X-ray binaries powered by accretion at super-Eddington mass-transfer rates (see Kaaret et al. 2017, for a review). They have long been thought to host intermediate-mass BHs (see e.g. Mezcua 2017), but the discovery of coherent pulsations from one of them (M 82 X-2; Bachetti et al. 2014) showed that ULXs can also be powered by accretion onto stellar-mass compact objects. Since then, five more pulsating ultraluminous X-ray sources (PULXs; Fürst et al. 2016; Israel et al. 2017a,b; Carpano et al. 2018; Rodríguez Castillo et al. 2020; Sathyaprakash et al. 2019) and a few candidates (e.g. Quintin et al. 2021) have been identified. At present, the ULX zoo is diverse, with the most luminous ones, called hyperluminous X-ray sources (Gao et al. 2003, L_X > 10⁴¹ erg s⁻¹), still consistent with accreting intermediate-mass BHs (e.g. ESO 243−49 HLX-1, Farrell et al. 2009) and those with luminosities in the range 10³⁹ − 10⁴¹ erg s⁻¹ consistent with accreting stellar-mass BHs or neutron stars (NSs), with the NS-ULX population possibly dominant (e.g. King & Lasota 2016; Koliopanos et al. 2017; Middleton & King 2017; Pintore et al. 2017; Walton et al. 2018a).

The spectral characteristics of ULXs suggest that observations may be influenced by two major factors: the mass-accretion rate (Sutton et al. 2013) and the inclination of the system with respect to the observer, with a possible precession of the accretion disc (Middleton et al. 2015a). According to the most accepted picture, matter accreting onto the compact object forms a disc with a complex structure: it starts as a standard thin accretion disc (Shakura & Sunyaev 1973) in the outermost part and becomes geometrically thick, inflated by the radiation pressure, until reaching its maximum height at the so-called spherisation radius, where the local radiation flux is equal to the Eddington flux. From the spherisation radius inwards, powerful outflows, which become clumpy at a certain distance, are launched from the surface of the disc (Shakura & Sunyaev 1973; Poutanen et al. 2007; Takeuchi et al. 2013). Following this idea, ULX X-ray spectra have been modelled with phenomenological components, with the aim to assess the different parts of the accretion disc. One possibility involves two multi-temperature disc models (Walton et al. 2018a,b; Gúrpide et al. 2021a, G21 hereafter), one with temperatures ≳1 keV, usually attributed to the innermost part of the accretion disc, and one with lower temperatures, attributed to the outermost regions of the disc (but see also Koliopanos et al. 2017, for a different interpretation of the components). Blueshifted absorption lines with velocities of a fraction of the speed of light at soft energies seen in a few ULXs (Middleton et al. 2014, 2015b; Walton et al. 2016a; Pinto et al. 2016, 2020, 2021; Kosec et al. 2018, 2021) would be the signature of the powerful ionised gas outflows. For now, the detection of these lines is at the limit of what current X-ray instrumentation is capable of, with detections at only the few sigma level.

On the basis of their spectral shape, ULXs are usually classified into three main categories: hard ultraluminous (HUL), soft ultraluminous (SUL), and broadened disc (BD; Sutton et al. 2013). The HUL and SUL sources can be modelled with a double black body, with temperatures around 0.1−0.3 keV and 1−3 keV for the cool and hot component, respectively. The dominance of one component with respect to the other determines the nature of HUL (hot component) or SUL (cool component) sources. On the other hand, the BD-ULX spectra can be modelled with a disc model with a radial temperature profile T(R)∝R^−p, with p < 0.75, indicative of an advection-dominated disc (Mineshige et al. 1994), as opposed to a standard thin disc, which would have p = 0.75 (Shakura & Sunyaev 1973). In addition, a sub-class of super-soft ultraluminous sources (SSUL) has also been identified, with ULXs exhibiting spectra dominated almost entirely by soft emission below ∼2 keV (see e.g. Urquhart & Soria 2016). Moreover, while some sources seem to permanently stay in the same class, others transit from one to another (Sutton et al. 2013; Middleton et al. 2015a; Gúrpide et al. 2021b).

M 81 X-6 is a bright ULX (L_X ≳ 2 × 10³⁹ erg s⁻¹) located in the outskirts of its host galaxy. Discovered first by Fabbiano (1988), this source was immediately classified as a ULX. The first observation of the source with the Chandra X-ray observatory (Weisskopf et al. 2000) was analysed by Swartz et al. (2003), who fitted the spectrum with a disc black-body model of ∼1 keV. Later, Gladstone et al. (2009) noticed a characteristic curvature above 2 keV and made the analogy of the spectral shape being similar to a BH binary in the hard state. Sutton et al. (2013) modelled the spectrum as a BD, but noticed that the source was still consistent with the HUL regime, especially considering its X-ray luminosity of ≳3 × 10³⁹ erg s⁻¹. Recently, Jithesh & Misra (2018) and Jithesh et al. (2020) conducted a broadband study of the source, analysing X-ray data from XMM-Newton (Jansen et al. 2001), NuSTAR (Harrison et al. 2013), and Suzaku (Mitsuda et al. 2007). Their broadband spectral analysis revealed similarities between M 81 X-6 and other PULXs, but a search for periodic modulation failed to identify any strong signal. They concluded that the source might be a potential PULX candidate. Finally, G21 explored this source in the context of a more general study on the long-term X-ray spectral evolution of ULXs. By tracking the evolution of the source in the hardness-intensity diagram (HID), they showed how the source goes back and forth between two main states. By comparing the evolution of the hardness with those of known PULXs, the authors identified M 81 X-6 as the best NS-ULX candidate of their sample. A similar bi-modal behaviour was also observed by Weng & Feng (2018) using observational data from the X-ray Telescope (XRT; Burrows et al. 2005) on board the Neil Gehrels Swift Observatory (Gehrels et al. 2004).

All of the aforementioned spectroscopic studies have used X-ray data taken up to 2012 and/or two simultaneous NuSTAR and Suzaku observations taken in 2015. Here, we aim to explore the evolution of M 81 X-6 with a more thorough and up-to-date analysis that accounts for all the available Chandra and XMM-Newton archival data from 2012 to 2020, complemented by all the data from Swift/XRT from its launch up to June 2021. In particular, this study benefits from several unpublished Chandra and XMM-Newton observations carried out in two different years (2016–2017 for Chandra and 2020 for XMM-Newton; see Table 1), which guarantee, for the first time, a spectral monitoring of the source over weekly, monthly, and yearly timescales. We also used the aforementioned NuSTAR dataset to investigate the spectrum at high energies (see Sect. 3).

2. Data reduction

2.1. Chandra

M 81 was observed by Chandra more than 15 times between February 2016 and January 2017. However, the source was within the field of view in only 11 observations, reported in Table 1. All observations were reprocessed with the Chandra Interactive Analysis of Observations (CIAO; version 4.13, Fruscione et al. 2006) software package, using the script chandra_repro, with calibration files CALDB version 4.9.5. For each observation, we extracted spectra and light curves, using circular extraction regions of 15″ for the source and of 45″ for the background, respectively. The unusually large radius for the source region was due to either the intensity of the source or to its off-axis position on the ACIS CCDs. We made sure that both the source and background extraction regions were on the same CCD to account for a correct subtraction of the background. Spectra were binned with a minimum of 20 counts/bin and light curves were extracted in three different energy bands: broad (0.3−10 keV), hard (1.5−10 keV), and soft (0.3−1.5 keV). Light curves were barycentred at the source coordinates RA = 9:55:32.9 and Dec = +69:00:33 (J2000, Gladstone et al. 2009) We also checked for pile-up, using the CIAO tool pileup_map¹, finding a level of pile-up < 5% in all the observations.

2.2. XMM-Newton

XMM-Newton archival data up to 2012 were already analysed in depth by G21. Therefore, we focused on more recent observations. XMM-Newton observed the galaxy M 81 14 times in 2020, but the ULX fell into the field of view in only eight observations (see Table 1). Six observations have no pn data because they were not centred on M 81 X-6 and the pn was in the small window or timing mode. XMM-Newton data were reduced with the Science Analysis System (SAS; version 19.1.0) and the respective current calibration files, and reprocessed with the standard SAS routines emproc/epproc, filtering the 10−12 keV light curve for count rates > 0.4 counts s⁻¹ to prevent contamination from the high-flaring background. Response files were generated using the tasks rmfgen and arfgen. We extracted spectra and barycentred light curves as described above for the Chandra data. Unfortunately, the XMM-Newton observations are much shorter than the Chandra observations, preventing us from using them for the timing analysis (cf. Sect. 4).

2.3. Swift/XRT

With the goal of exploring the gap between the latest Chandra and XMM-Newton observations and those in G21, we looked into the Swift data archive². For each of the 603 observations, we used the online tools³ (Evans et al. 2007, 2009) to build the light curves in the same broad, hard, and soft energy bands to retrieve the hardness of the source. Figure 1 shows the broadband (0.3−10. keV) Swift/XRT light curve, with the times of Chandra, XMM-Newton, and NuSTAR observations highlighted.

Fig. 1. Temporal evolution of M 81 X-6. Left panel: Swift/XRT light curve in the 0.3−10 keV energy band. Each data point corresponds to an observation. The time interval spun by the Chandra observations considered in this work is coloured in light blue, while that of the most recent XMM-Newton observations is coloured in light orange or represented by a dotted orange vertical line. The vertical dotted green lines mark the oldest XMM-Newton observations (cf. G21) and the red one the NuSTAR observation used in this work. Right panel: histogram of Swift/XRT count rates, showing a clear bimodal distribution.

2.4. NuSTAR

NuSTAR has observed the galaxy M 81 in seven instances, catching the ULX only in one of the observations (ObsID 60101049002). The data were extracted with the NuSTAR Data Analysis Software version 2.0.0 with CALDB version 20200826. Source and background spectra were extracted using the task nuproducts with the standard filters. Source events were selected from a circular region centred on the source of 60″ in radius, to avoid contamination from the nearby active galactic nucleus in M 81. Background regions were selected from larger, circular, source-free regions and on the same chip as the source but as far away as possible to avoid contamination from the source itself. NuSTAR spectra were re-binned to 40 counts per bin due to its lower energy resolution and the lower signal-to-noise ratio compared to the EPIC cameras.

3. Spectral analysis

As mentioned in Sect. 1, M 81 X-6 shows spectral changes between states at high and low luminosity. To better understand and eventually classify the source spectral states, as a first step, we decided to follow the method of G21 and fit each individual Chandra and XMM-Newton spectrum with two double-absorbed disc black bodies. All fits were performed with the X-ray spectral fitting package Xspec (Arnaud 1996), version 12.12.0b, with abundances from Wilms et al. (2000) and cross-sections from Verner et al. (1996). The two black-body components were reproduced by the Xspec model diskbb (Makishima et al. 1986; Mitsuda et al. 1984), so that the overall fit function was: tbabs*tbabs*(diskbb+diskbb). The Galactic absorption was fixed to 9.8 × 10²⁰ cm⁻² (HI4PI Collaboration 2016) and local absorption was fixed to 11.7 × 10²⁰ cm⁻², the average value determined by G21, since it was not always possible to constrain it. XMM-Newton observations where data from more than one instrument (MOS1, MOS2 and/or pn) were present simultaneously were fitted together. We did not perform a spectral fit only for the XMM-Newton ObsID 0870930401, due to the low number of counts. Best-fit values are reported in Tables 2 and 3, with errors at 90% confidence level.

As already observed in G21 for M 81 X-6, as well as for several other ULXs, the fit with two thermal black-body components is in general a good fit, with reduced χ² values close to 1, except for a few observations, where it is lower than 1 (noticeably Chandra ObsIDs 18052 and 19993, and XMM-Newton ObsIDs 0870930101, 0870930501, and 0870931001, most likely because of the low number of counts). Unfortunately, because of the quality of the data, for several Chandra observations the best-fit values were not always constrained and in many cases we were only able to place lower limits. Overall, the best-fit temperatures of the soft component are around 0.4 keV, but while they are well constrained for XMM-Newton, those of Chandra have very large uncertainties. The best-fit values of the hard component are always higher than 1.3 keV. Both best-fit temperatures are in general consistent with those found by G21 for the 2001−2012 XMM-Newton datasets, as well as the unabsorbed luminosities in the range ∼2 − 5 × 10³⁹ erg s⁻¹.

We noticed that in almost all Chandra observations and in two XMM-Newton ones (ObsID 0870900101 and 0870900201) residuals at energies ≲1 keV are marginally evident. Those residuals have sometimes been interpreted as signatures of powerful outflows irradiated from the surface of the accretion disc (Middleton et al. 2014, 2015b). In order to resolve the spectroscopic features of these winds, high-resolution spectra are necessary (Pinto et al. 2016, 2020, 2021; Kosec et al. 2021). Hence, we did not add any additional feature to the fit model.

As in G21, we also computed the hardness ratio (HR) as the ratio between the unabsorbed flux in the hard band (1.5−10 keV) over the unabsorbed flux in the soft band (0.3−1.5 keV) and plotted it in the HID of Fig. 2, where we also show the results from G21 for an easy comparison. Our results are consistent with those of G21, except for a single old XMM-Newton observation, taken in 2004, where the source was caught with a higher luminosity (∼8 × 10³⁹ erg s⁻¹), while in the recent observations it never went above L_X ∼ 5 × 10³⁹ erg s⁻¹.

Fig. 2. Hardness-intensity diagram (HID) of the Chandra (blue) and 2020 XMM-Newton (orange) datasets, in comparison with the old XMM-Newton observations (fainter text and points) as reported in G21.

We also investigated whether M 81 X-6 entered other regions of the HID, as seen in Holmberg II X-1 and NGC 5204 X-1 (Gúrpide et al. 2021b). To this end, we considered all Swift/XRT archival data up to June 2021. Following Gúrpide et al. (2021b), we computed the HR for each observation as the ratio of the counts in the hard (1.5−10 keV) and soft (0.3−1.5 keV) energy bands. To directly compare Swift/XRT data with those from the other two X-ray satellites, we converted Chandra and XMM-Newton count rates into Swift/XRT count rates. Briefly, we convolved the absorbed best-fit models of each Chandra/XMM-Newton observation with the Swift/XRT redistribution matrix (swxpc0to12s6_20130101v014.rmf) and ancillary file (swxs6_20010101v001.arf) of 01/08/2020. Results are shown in Fig. 3, where the broadband (0.3−10 keV) count rates are plotted against the HR. Chandra and XMM-Newton data are perfectly consistent with Swift/XRT with no overall evidence of a third state, where the source is supposed to be harder and/or at higher luminosities. Our analysis confirms the results of Weng & Feng (2018).

Fig. 3. Evolution of Swift/XRT broadband count rates (blue dots) as a function of the HR, in comparison with the latest Chandra (empty orange circles) and XMM-Newton data (empty orange diamonds). Data from different XMM-Newton-EPIC instruments are plotted individually. Errors on Chandra and XMM-Newton data are not shown.

From the inspection of the individual Chandra and XMM-Newton spectra and of the HIDs, we found that the source transits between two main regimes: a hard, high-luminosity (HHL) state, characterised by an HR ∼ 2−3 and luminosities of the order of 4 − 6 × 10³⁹ erg s⁻¹, and a soft, low-luminosity (SLL) state, with HR ≲ 2 and L_X ∼ 3 × 10³⁹ erg s⁻¹ (see Tables 2 and 3). Figure 4 shows two spectra representative of each group. We note that these two regimes do not fit in the traditional HUL and SUL states (cf. Sect. 1). The overall hardness of the source is always ≳ 1 (see e.g. Figs. 2, 3, and 5), which would classify M 81 X-6 as a HUL source. We simply note that it is within the same HUL category that M 81 X-6 explores two different regimes. These two regimes are also evident in the double-peaked distribution of Swift/XRT count rates (Fig. 1, right panel).

Fig. 4. Spectra of Chandra ObsID 18047 (blue) and XMM-Newton/pn ObsID 0870900201 (orange), representative of the HHL and SLL states, respectively. Spectra are fitted with the model tbabs*tbabs*(diskbb+diskbb) (solid line) with the two black-body components plotted as dashed lines. Bottom panel: residuals of the fit.

Fig. 5. Light curves and HR of Chandra ObsID 18047. Top panel: light curves in the hard (orange) and soft (blue) energy bands, binned at 1000 s, with dashed lines representing the mean count rate of each light curve. Bottom panel: HR, computed as the number of counts in the hard band divided by the total number of counts. The dashed line corresponds to the mean HR.

Once we had established the boundaries of the two spectral states seen for M 81 X-6, we grouped together the spectra of the same type and fitted them, with the aim to constrain the temperatures and the normalisations of the two black bodies. At this stage, to obtain statistically acceptable fits, we had to exclude Chandra ObsIDs 18875 (HHL) and 18054 (SLL) and XMM-Newton ObsIDs 0870900101 and 0870930401 (SLL), because they represented transitory stages between the two regimes. Best-fit results of the state-resolved spectra are reported in Table 4. The best-fit temperatures resulted in keV and keV for the HHL regime, and keV and keV for the SLL regime. These results would suggest that the soft component does not vary between the two regimes, in agreement with G21. Thanks to the high number of total counts, this time we were able to constrain also the local absorption, which we left free to vary. It resulted in cm⁻² and cm⁻², for the HHL and the SLL regimes, respectively, consistent in both cases with the average value of 11.7 × 10²⁰ cm⁻² estimated by G21.

Whilst these last fits are statistically acceptable, we noticed that the best-fit temperature of the hard component is higher in the SLL state rather than in the HHL state. This may imply that the hard part of the spectrum is not well-described by the dual thermal component model, as suggested by weak residuals at energies ≳8 keV seen in a few spectra (e.g. ObsIDs 18054, 18047, 18817, 19685). In the literature, an additional cut-off power law has often been used by several authors to model high-energy excesses and it has been linked to the X-ray radiation from the accretion column of the NS in PULXs (e.g. Walton et al. 2018a, cf. Sect. 5). Hence, we decided to add this further cut-off power law component (cutoffpl in Xspec) to the best-fit model, thus allowing for the possibility of a NS accretor also for M 81 X-6. To constrain the parameters of the cut-off power law, we used the NuSTAR observation, taken when M 81 X-6 was in the SLL regime (see Fig. 1). Hence, at first, we fitted all the XMM-Newton observation in the SLL state used in the previous step together with the NuSTAR/FPMA and FPMB data, with an absorbed double disc black-body plus cut-off power law model (i.e. tbabs*tbabs*(diskbb+diskbb+cutoffpl)). NuSTAR data were limited to the 3−20 keV range, because of the high background contamination at higher energies. To compute the errors on E_cut and Γ, we froze the local absorption and alternately one of the two parameters, leaving the other free to vary. We found the best-fit values of the high-energy cut-off to be keV and for the power law index . These values are consistent with those of other ULXs (see e.g. Walton et al. 2018a). To assess whether adding a cut-off power law actually improves the fit, we performed an F-test, which resulted in a statistical value of 4.2 and a probability of 0.04, meaning that the additional cut-off power law component is marginally significant. Next, we froze the power law parameters to their best-fit values and fitted all the spectra in the two regimes with the same model, with the local absorptions, the two black-body temperatures and all the normalisations tied for each group. Best-fit results are shown in

Table 4. We also performed an F-test for the HHL dataset, which was inconclusive. This is not surprising, since the cut-off power law is well characterised for the SLL state but not for HHL, due to the lack of high-energy coverage of this latter state. Nonetheless, we chose to keep the cut-off power law also for the HHL state, in light of the interpretation of it modelling the X-ray emission coming from the accretion column of the NS. As a matter of fact, because we expect the cut-off power law to be stronger in the HHL state, the best-fit values found for the SLL state could be considered lower limits for the parameters of the cut-off power law in the HHL state. In any case, whether we take the cut-off power law into account or not, the best-fit values of the black-body components of the HHL state remain unchanged. Overall, we found the best-fit temperatures of the soft black-body component of the two regimes in agreement with each other (T_soft ≈ 0.4 keV), while the best-fit temperature of the hard black body seems to be higher in the HHL state ( keV) and lower in the SLL state ( keV), though they are consistent within the uncertainties.

To further constrain the soft component, we performed one last, simultaneous fit of the two groups, where we linked the temperature and the normalisation of the soft component for all the observations. We left free to vary within each group the local absorption, the temperature and normalisation of the hard component, and the normalisation of the cut-off power law. We obtained a fit with (d.o.f.) = 1.06(1985) and a null hypothesis probability (n.h.p.) of 0.03. The best-fit temperature of the soft component was keV, with a normalisation of . As expected, this best-fit temperature is consistent with the best-fit temperatures of the soft component of the two states fitted separately. All the other best-fit parameters were also consistent with the previous fits. Although the value of the n.h.p. would imply that this fit is only marginally statistically acceptable, the consistency of all the best-fit parameters with the previous fits reassures on its validity. Nevertheless, in the successive steps of our analysis we consistently use the best-fit of T_soft ∼ 0.4 keV, as obtained in the separate fits of the two regimes.

We also studied the relationship of the hard and soft temperatures with their unabsorbed bolometric luminosities (in the range 0.01−100 keV). Results are shown in Fig. 6, together with values from G21. We did not find any trend, neither for the hard nor the soft component. Nonetheless, we plot the expected L − T relations of L ∝ T⁴ in the case of a thin disc (Shakura & Sunyaev 1973) and L ∝ T² for an advection-dominated disc around a BH (see e.g. Watarai et al. 2000; Walton et al. 2020), as well as the trend L ∝ T⁻⁴ recently proposed for NGC 1313 X-2 (Robba et al. 2021). The lack of any trend for our dataset is not surprising, as we are attempting to use simple models to fit ULX spectra when the physics involved is certainly more complex. Moreover, as remarked in G21, these correlations are likely due to the degeneracy between the temperature and the normalisation of each component, casting doubts on their physical interpretations.

Fig. 6. Bolometric luminosity as a function of the temperature for the soft (left panel) and hard (right panel) black-body components, in comparison with the older XMM-Newton observations as in G21 (in grey). For the soft component, three L − T laws are overplotted: L ∝ T4 for a thin disc (Shakura & Sunyaev 1973), L ∝ T2 for an advection-dominated disc around a BH (Watarai et al. 2000), and L ∝ T−4 as observed for NGC 1313 X-2 (Robba et al. 2021). The bolometric luminosities are computed in the 0.01−100 keV range.

As mentioned in Sect. 1, the HHL state of M 81 X-6 has sometimes been also modelled with a BD component (e.g. Sutton et al. 2013), such as diskpbb in Xspec (Mineshige et al. 1994). To test this model, we tried to substitute the harder multi-temperature disc component with a BD for the spectra in HHL state, so that the overall model is now: tbabs*tbabs*(diskbb+diskpbb+cutoffpl). We note that the diskpbb alone can likely reproduce the entire 0.3−10 keV spectrum, but, because we wanted to test a model with two thermal components, we kept them both and froze the parameters of the soft component to ensure its significant contribution to the fit. In this way, we obtained a best-fit temperature of the BD component of keV, a parameter , and a normalisation of , with a (d.o.f.) = 1.06(1151) and a n.h.p. of 0.08. The parameter p is consistent at 90% confidence level with the value of 0.75 for a thin Shakura-Sunyaev disc. Hence, there is only marginal evidence of the hard disc component to be advection-dominated. We found similar results for the SLL spectra. For this reason, we do not investigate further and consider the model with two disc black bodies plus the cut-off power law as the best-fit model.

4. Temporal analysis

As for other ULXs, the X-ray light curves of M 81 X-6 show short-term variability ranging from a few tens to thousands of seconds. Such variability is stronger in the hard band than in the soft one, as seen in many ULXs (e.g. Middleton et al. 2011; Sutton et al. 2013; Weng & Feng 2018). An example is shown in Fig. 5, which also displays the HR as the fraction of counts in the hard band over the total counts in the broad energy band, for an easier interpretation. The same trend was also observed in Swift/XRT data by Weng & Feng (2018) (cf. Sect. 5) and it is consistent with the lack of variability of the soft spectral component presented above.

To quantify the overall degree of variability of the source, we computed the fractional variability F_var (i.e. the square root of the normalised excess variance) as in Vaughan et al. (2003) for the Chandra data, dividing the light curves into segments of 20 ks, corresponding to the exposure of the shortest observation, and averaging the obtained values. Due to the good time intervals and the intrinsic quality of the data, only five datasets allowed for a correct estimation of the fractional variability (see Table 5). All five values are consistent with one another within the errors, so no trend of F_var with hardness or luminosity can be established (see also Fig. 7).

Fig. 7. HID of the Chandra observations and one XMM-Newton observation taken in 2001 (G21). The circle size is proportional to the fractional variability, only shown on observations for which it could be constrained.

As stated in Sect. 2.2, XMM-Newton observations are too short for a meaningful timing analysis. However, we did compute the fractional variability for the old XMM-Newton ObsID. 0111800101, taken in 2001 (see G21), which has a good exposure of 73.6 ks for the pn. It resulted in F_var = 0.037 ± 0.005, computed over four 20-ks intervals, and hence consistent with the values of the Chandra observations.

Unfortunately, with the present data, we could not conduct a search for pulsation of the order of seconds, as found for PULXs, since none of the XMM-Newton/pn observations has a sufficient number of counts (see Rodríguez Castillo et al. 2020). Nonetheless, we determined the upper limit on the pulsed fraction (PF) for the XMM-Newton/pn observation of 2001. We used the software HENDRICS, version 7.0.2 (Bachetti 2018) on the barycentred events at the Chandra position of M 81 X-6 (Swartz et al. 2003) on the 0.3−12 keV band and searched for purely sinusoidal profiles with the statistic (i.e. n = 1), using 12 bins for each folding frequency and scanning the 0.1−8 Hz frequency range (similar to other ULXs). We used the option fast, which automatically optimises the f − ḟ search space. We did not detect any signal in any of the long ∼20 ks good time intervals (GTIs), and placed an upper limit on the PF ∼10.5% at a 90% confidence level. Combining three of the long continuous GTIs yielded a more stringent upper limit of ∼7.2%. We did not compute PF upper limits for the other XMM-Newton/pn observations, due to their short exposures.

5. Discussion

We have presented a thorough analysis of all available X-ray data of the ULX M 81 X-6 taken with Chandra, XMM-Newton, Swift/XRT, and NuSTAR in the last decade. Older XMM-Newton archival observations, up to 2012, were already analysed by G21. Here, we built on their work and presented an up-to-date analysis of the source.

5.1. M 81 X-6 spectral states

According to the general classification of ULXs (Sutton et al. 2013), with its HR always ≳1 (see Figs. 2 and 3), M 81 X-6 would fall into the general category of HUL sources. However, as already noted in G21, within its HUL state it oscillates between two main regimes, which we defined as HHL (HR ∼ 2−3 and L_X ∼ 4 − 6 × 10³⁹ erg s⁻¹) and SLL states (HR ≲ 2 and L_X ∼ 3 × 10³⁹ erg s⁻¹; see Fig. 2 and Sect. 3). A comparison with the much more extended Swift/XRT data sample (Fig. 3) confirms that the source oscillates over the years between these two main states, without any hint of any further state, as opposed to, for example, Holmberg II X-1 and NGC 5204 X-1, (Gúrpide et al. 2021b). The two different types of track found for the ULXs in the sample of G21 led the authors to suggest that the nature of the accretor (BH vs NS) may be responsible for the spectral differences. As a matter of fact, the track on the HID of M 81 X-6 is similar to those of other ULXs (e.g. M 83 ULX1, Circinus ULX5, Holmberg IX X-1, NGC 55 ULX1, and M 51 ULX8; see G21) and remarkably similar to the PULX NGC 1313 X-2 track (see also Sathyaprakash et al. 2022). It also shares the same range of hardness (HR ∼ 1−3) and luminosity (L_X ∼ 2 − 10 × 10³⁹ erg s⁻¹) with NGC 1313 X-2, leading to the idea that M 81 X-6 could host a NS as well. Under the hypothesis of a NS central engine, a lack of detected pulsation can be attributed to the relative inclination of the rotational axis of the system and the direction of the line of sight (G21), to the scattering within the wind cone (e.g. Mushtukov et al. 2021), or to the relative position of the magnetospheric and the spherisation radii (e.g. Walton et al. 2018a). A separate analysis of the two states would hence shed light on the geometry of the system and on the physical processes at play, as we discuss below and in the following sections.

We sorted all the Chandra and XMM-Newton observations by state and fitted the two groups, to better constrain the parameters of the model. The simple model of two thermal discs that we tried at first resulted in a scenario in which the temperature of the hard component was higher in the SLL state than in the HHL state. Motivated by this possible contradiction and by the presence of small residuals at energies ≳8 keV, we decided to add a further component to the fits, to better model the harder part of the spectra. To this purpose, we chose a cut-off power law, which we constrained thanks to the NuSTAR observation, taken when the source was in the SLL regime. With this new model, we found a temperature of the soft disc component equal to ∼0.4 keV in both states, while the temperatures of the hard disc component are ∼1.5 keV and ∼1 keV, in the HHL and SLL states, respectively. We must stress once again that the lack of any NuSTAR observation of M 81 X-6 in the HHL at high energies prevented us from correctly modelling the cut-off power law in this state. Nevertheless, we resolved to use this component to fit the HHL spectra as well, under the assumption that the disc temperatures are more meaningful.

5.2. M 81 X-6 as a NS-ULX

Walton et al. (2018b) found that the spectra of the PULXs NGC 5907 ULX and NGC 7793 P13 show a hard excess above ∼10 keV when fitted with two thermal disc models. Based on phase-resolved spectroscopy, the authors determined that this hard excess was attributed to the pulsed component, modelled with a cut-off power law. Moreover, Walton et al. (2018a) found that those PULXs had a larger contribution to the 0.3−40 keV flux from this component (typically above 50%), compared to ULXs for which pulsations had not been found. In light of that, we computed the ratio between the flux of the cut-off power law and the total flux F_cutoff/F_tot (Table 4), using the same energy range of Walton et al. (2018a). We found a ratio of 50% in the SLL regime, where we had NuSTAR coverage. It is interesting to note that the value found for the SLL sits in between those seen in PULXs and the rest of the ULX population with broadband spectroscopy available. Given the similarities of M 81 X-6 with NGC 1313 X-2 and the low PF of the latter (the lowest in the PULX population), our analysis would support that M 81 X-6 indeed hosts a NS with a small PF. We also found that the F_cutoff/F_tot ratio drops to below 2% in the HHL regime. This result is more uncertain, given that we had to assume the same parameters for the cut-off power law as for the SLL state, due to the lack of high-energy (> 10 keV) coverage in this state. Nonetheless, it suggests that pulsations should be more readily observable in the SLL state, similarly to NGC 1313 X-2, where pulsations were detected in the low-luminosity state (Sathyaprakash et al. 2019, G21). Moreover, the PF for the XMM-Newton observation taken in 2001 resulted in an upper limit of ≲10% (cf. Sect. 4), consistent with the PF of other ULXs (Fürst et al. 2016; Rodríguez Castillo et al. 2020), but higher, for instance, than the PF of NGC 1313 X-2 (Sathyaprakash et al. 2019). Whether this low PF is specific to this spectral state or applies to all states remains to be seen.

5.3. Long-term variability

M 81 X-6 shows a small variability of the broad-band light curves throughout its evolution, with low values of fractional variability (F_var ≲ 10%; see Table 5 and Fig. 7), close to those of NGC 1313 X-2 (Robba et al. 2021), for which a correlation of the fractional variability with the luminosity has been observed (it has to be noted, however, that other ULXs show higher values; see e.g. Middleton et al. 2015a; Salvaggio et al. 2022). Unfortunately, in our case the quality of the data and/or the short length of the observations, especially those made with Chandra, meant large error bars for the fractional variability so that it is difficult to infer a trend within the two states. Nonetheless, one of the oldest XMM-Newton/pn observations (ObsID. 0111800101; G21) was long enough to constrain its fractional variability with a sufficiently high precision. This particular observation took place when the source was in a transitional state close to the HHL regime, but its low fractional variability of 0.037±0.005 alone would seem to contradict the trend found for NGC 1313 X-2. More observations would be of interest to confirm or rule out this trend.

Nevertheless, M 81 X-6 exhibits a bi-modal distribution on the HID, already noticed by Weng & Feng (2018) for Swift/XRT data up to August 2017 (cf. Fig. 1, right panel), with a strong quasi-periodic modulation in the hard band (1−10 keV), but not in the soft one (0.3−1 keV). These spectral changes seem to be associated with a quasi-periodicity of ∼110 days (Lin et al. 2015). Lin et al. (2015) also reports another periodicity at ∼370 days, but we note it is uncertain as it was only supported by two cycles and, given the slow timescale, could be a spurious feature. This bi-modal behaviour, common to other investigated ULXs (e.g. Holmberg IX X-1 and X-2), has been attributed to precession of the accretion flow, responsible for the modulation of the hard X-rays coming from the innermost regions, while the isotropic soft emission arising from the outflows and the photosphere of the outer disc would remain more or less constant (cf. Sect. 5.4).

Alternatively, a bi-modal distribution, with variations of more than an order of magnitude in luminosity observed for other ULXs, has been attributed to the propeller effect (Earnshaw et al. 2018; Song et al. 2020). However, the difference in luminosity of M 81 X-6 between the HHL and the SLL regimes is slightly less than a factor of two. The ratio of the luminosities in/out of the propeller regime can be estimated as (Tsygankov et al. 2016, Eq. (5))⁴, where P is the NS rotational period, M_1.4 its mass in units of 1.4 M_⊙, and R₆ it radius in units of 10 km. Using the canonical values for a NS of 1.4 M_⊙ and 10-km radius, a quick estimation would lead to P ∼ 1.28 ms, far too small compared to the spin period of known PULXs. Hence, assuming that all NS-ULXs share the same properties as those already observed, we can exclude propeller transitions as the mechanism responsible for the variability. This conclusion is obviously based on a rather small sample of PULXs and on the limitations of current X-ray instruments in detecting those with shorter periods. Increasing the sample of PULXs will be crucial to properly investigate the role of the propeller regime in ULXs.

5.4. The nature of the super-orbital period

If we consider a scenario in which the soft thermal component represents the outflow photosphere, where the temperature is a function of the mass-accretion rate (Poutanen et al. 2007), then the lack of variation in T_soft between the two regimes suggests that the luminosity and spectral variations of the source are not driven by changes in the mass-transfer rates, in line with G21. The bi-modal variability and the presence of a variable ∼110 days quasi-periodicity (Lin et al. 2015) may therefore imply that these changes are driven by precession of the accretion disc. The nature of this super-orbital period in ULXs is still uncertain (e.g. Vasilopoulos et al. 2020). Here, we explore two possible scenarios: Lense-Thirring precession and the torque from the NS magnetic dipole. We did not explore free precession (see e.g. Vasilopoulos et al. 2020), because the variability in the super-orbital periods of some sources (e.g. M 51 ULX7; Brightman et al. 2022) would seem to disfavour this scenario.

Lense-Thirring precession is caused by the misalignment of the rotational axis of the system and NS spin axis, which induces frame-dragging and oscillations in the disc. Middleton et al. (2019) suggested that BH- and NS-ULXs should be distinguishable on the basis of the temperature from the spherisation radius (R_sph) and super-orbital period, given the physical differences between these two objects (e.g. spin, mass, presence of a surface, etc.). Additionally, it may be possible to offer some constraints on the magnetic field of the source (B) and its spin (P_NS). To this end, under the assumption that Lense-Thirring precession drives the variability of the source, we employed the formulae presented by Middleton et al. (2019; see also e.g. Poutanen et al. 2007; Mushtukov et al. 2019; Vasilopoulos et al. 2019; Middleton et al. 2018; Fragile et al. 2007), setting all relevant unknown parameters (such as the accretion efficiency or the ratio of the asymptotic wind velocity relative to the Keplerian velocity, β) as in the original work of Middleton et al. (2019), unless stated otherwise. We note that we assume the precession to extend until the photospheric radius (as for Middleton et al. 2019), but this should be regarded as an upper limit for the radius at which the flow is able to precess as a solid body, as the bending wave might not be able to propagate if the outflow becomes supersonic at smaller radii (see Middleton et al. 2019).

A complete modelling of the Lense-Thirring precession requires the tuning of two more parameters: the dimensionless NS spin, a_*, and the fraction of energy that powers the outflow, ε_wind = L_wind/(L_wind + L_rad), with L_wind the kinetic luminosity of the wind and L_rad the observed radiative luminosity (Poutanen et al. 2007). Based on the periods found in ULXs (e.g. Bachetti et al. 2014; Israel et al. 2017a; Rodríguez Castillo et al. 2020), a_* should be ≪0.01, as already pointed out by Middleton et al. (2019). Assuming canonical values for the moment of inertia of the NS (I = 10⁴⁵ g cm²) and its mass (M = 1.4 M_⊙), we find that 3 × 10⁻² > a_* > 3 × 10⁻⁵ should approximately cover the range of observed periods in ULXs (∼0.1 s < P < 10 s). For ε_wind, we considered values of 0.2 and 0.8 to cover the range of values derived from numerical simulations (e.g. Jiang et al. 2014) and some observational constraints derived from shock-ionised optical bubbles around ULXs, as for, for example, NGC 1313 X-2 (Pakull & Mirioni 2002), IC 342 X-1 (Cseh et al. 2012), NGC 5585 X-1 (Soria et al. 2021), and NGC 1313 X-1 (Gúrpide et al. 2022).

The resulting Lense-Thirring plane is shown in Fig. 8, together with the position of M 81 X-6 in it and the expected values for a 20 M_⊙ BH with a moderate spin of 0.1 for reference. We note that attributing the temperature derived from our simple modelling to the temperature from R_sph is certainly problematic as (a) our modelling is not able to capture the complexity of the accretion flow and (b) it is unlikely that the hottest temperature returned by the diskbb from a range of black bodies coincides with the temperature from R_sph. Because of the latter reason, our naive expectation would be that our estimate for T_sph constitutes an upper limit. These caveats should be borne in mind when interpreting these results and therefore our estimates should be consider semi-quantitative. Improvement in theoretical modelling of PULXs should help reduce these uncertainties in the future.

Fig. 8. Lense-Thirring accretion plane adapted from Middleton et al. (2019). Left and right panels: resulting Tsph and precessing period of the wind for εwind = 0.2 and 0.8, respectively. The values expected for a 20 M⊙ BH with a dimensionless spin of 0.1 are shown in dashed dark blue and could be considered as a rough upper limit for the location of BHs in the plane. The rest of the lines show the values for NSs with different dipole magnetic field strengths, as indicated. The solid and dashed lines show the effect of varying the spin of the NS from 3 × 10−5 to 3 × 10−3, which we have converted to periods assuming canonical values for the moment of inertia and mass of the NS (see the text for details) and represent typical values observed in ULXs.

From Fig. 8, we see that the position of M 81 X-6 on the plane indicates a weakly magnetised NS (B < 10¹⁰ G) as the accretor. We note that lower magnetic field values would yield equivalent results for M 81 X-6 as those found for B = 10¹⁰ G. This is because for low magnetic field values (B < 10¹⁰ G) and accretion rates with T_sph ≲ 0.6 keV we find already that (i.e. the magnetospheric radius R_mag is pushed all the way to the innermost stable orbit R_ISCO) and therefore the solutions become independent of the magnetic field. This implies that our magnetic field estimate constitutes an upper limit, which could be as high as 10¹² G if we consider the curves for ε_w = 0.2. We note however that mostly B < 10¹⁰ G would seem favoured as ε_w = 0.2 constitutes a limiting case and high values would rule out the case for B ∼ 10¹² G. Moreover, since higher values of ε_w have the effect of increasing the size of the magnetosphere, as more mass is lost into the outflow, for ε_w > 0.3, a magnetic field with B = 10¹² G in M 81 X-6 would truncate the disc before becoming supercritical (this effect can be seen on the Fig. 8 as the regions where the lines cut, in particular on the right panel), which would be inconsistent with this scenario.

Therefore, according to this scenario and considering the super-orbital period P = 111 d (Lin et al. 2015), T_sph ∼ 0.4 keV, ε_wind = 0.2 and using B = 10¹⁰ G, we would expect for M 81 X-6 a spin period P_NS ∼ 4.5 s, a mass-accretion rate at the magnetospheric radius ṁ₀ = 42, a mass-accretion rate at the inner disc radius ṁ_in ∼ 37, a magnetospheric radius R_mag ∼ 13 R_g, a spherisation radius R_sph ∼ 302 R_g, and an outer atmosphere radius R_out ∼ 758 R_g, where R_g is the gravitational radius (R_g = GM_NS/c²). For ε_wind = 0.8, we found P_NS ∼ 4.7 s, ṁ₀ = 21, ṁ_in ∼ 12, R_mag ∼ 18 R_g, R_sph ∼ 129 R_g, and R_out ∼ 726 R_g. The first thing to note is that the predicted spin periods for the NS are within the range of those observed in ULXs and therefore not implausible. The values of ṁ₀ and B found for both ε_w seem to be within the allowed parameter space for the Lense-Thirring torque to dominate over the tidal torque according to Middleton et al. (2019, their Fig. 3). Therefore, it might be reasonable to suggest that this torque might be at play in this source, and even be responsible for the precession. Finally, this scenario would support the identification of M 81 X-6 as a weakly magnetised NS, where R_sph > R_mag, as suggested by G21. The overall low R_mag, close to the R_ISCO = 6 R_g, coupled with the large R_sph might explain the lack of detectability of pulsations, due to a low PF component as a result of the small extent of the magnetised region or to the signal being scattered and lost in an accretion curtain (Mushtukov et al. 2017).

Given the relatively slow P_NS required for this scenario, however, it is possible that the magnetic torque due to the interaction from the NS magnetic field with the accretion disc (Lipunov & Shakura 1980) will dominate over the Lense-Thirring torque, notably for magnetic fields ≳10¹¹ G (see Middleton et al. 2018, their Fig. 8). The exact calculation is not straightforward as the magnetic torque has only been estimated for a standard thin disc (see Middleton et al. 2018, and references within), but it is informative to consider whether this scenario could account for the super-orbital period of the source, as also explored by Vasilopoulos et al. (2020) in the case of the PULX M 51 ULX7 (see also Mushtukov et al. 2017). The precessing timescale of the NS can be estimated as (Lipunov & Shakura 1980)

(1)

where ψ is the angle between the NS rotational axis and the normal to the accretion disc and ξ is the angle between the NS magnetic field and its rotational axis. R_NS is the radius of the NS at its surface, assumed to be 10 km. As above, we solve numerically for R_mag, based on ṁ₀ estimated from the temperature of the soft component and for different NS magnetic fields. We set ξ = ψ = 30° for reference and carry out the calculations for ϵ_w = 0.2 and 0.8 as above and for P_NS = 0.1 − 20 s. We present the calculations in Fig. 9, which can also be useful for future reference.

Fig. 9. Precessing period due to the interaction of the NS magnetic field with the accretion disc (Lipunov & Shakura 1980) as a function of the NS spin period calculated for different values of the NS dipole magnetic field. Left and right panels: calculations for ϵw = 0.2 and ϵw = 0.8, respectively. The horizontal dashed line indicates the super-orbital period of 110 days of M 81 X-6 (Lin et al. 2015), and we have considered the temperature of the soft component in the calculations of Rmag.

As opposed to the Lense-Thirring scenario, this model would require a strongly magnetised (10¹² G < B < 10¹⁴ G) NS for M 81 X-6, whereas the necessary spin period would be in line with those estimated under the assumption of Lense-Thirring precession. We thus see that estimating the magnetic field of the putative NS could be useful in putting constraints on the mechanism responsible for the super-orbital period.

It must be noted that if this mechanism operates in the whole PULX population, it would require magnetic fields well above 10¹⁵ G for those more rapidly spinning PULXs and/or with shorter super-orbital periods (e.g. M 51 ULX7; as already pointed out by Vasilopoulos et al. 2020). Therefore, if the same mechanism is not responsible for the super-orbital periods seen in other PULXs (e.g. Walton et al. 2016b), then we may be able to rule out this scenario for M 81 X-6 too. We conclude that, from the mechanisms explored here, Lense-Thirring precession can explain the super-orbital period of M 81 X-6, provided that the NS has a magnetic field ≲10¹⁰ G.

Finally, to come full circle, the characteristic radii derived above can be compared with the ‘apparent’ inner emitting radii (Kubota et al. 1998) of the two multi-temperature disc components. These can be computed from the normalisations of the two black bodies, for which we assume two extreme inclination angles of 0° and 30° inferred in G21, as the source is viewed at low inclinations (see also G21). Results are shown in Fig. 10. In our interpretation of the system, the soft emission is arising from the outer disc and/or the outflows, between R_sph and R_out, whereas the harder emission would come from the inner region from within R_sph down to R_ISCO for a BH-ULX or to R_mag for a NS-ULX. The temperature of the soft disc component, which corresponds to the temperature of the disc at the spherisation radius, has been directly employed in the Lense-Thirring model. Then, it is not surprising that the apparent inner radii obtained from the soft component are consistent with R_sph or fall within R_sph and R_out. The only exception may be the case in the HHL state with ε_wind = 0.2 (Fig. 10, top left panel), where the apparent inner radius falls slightly before R_sph, but, in any case, the soft-emission region would still develop in between the R_sph and the R_out. On the other hand, we find relatively good agreement between R_mag and the radii derived from the hard emitting component, even though this component is not directly implied in the Lense-Thirring model. In all possible scenarios, the inner radii fall within R_mag and R_sph. This tentative agreement may suggest that our simplistic spectral fitting model might be a good first order approximation and may strengthen our interpretation of the various emitting regions attributed to the spectral components. A more accurate comparison should take into account, for example, the fraction of energy deposited in the wind, the colour-temperature correction, the fact the photosphere is anistropic, etc., but this advanced modelling is far beyond the goal of this paper.

Fig. 10. Comparison of the physical sizes of the system derived from the best-fit multi-temperature disc components (dots) and from the Lense-Thirring precession model (vertical lines; see Sect. 5.4 for details). The solid and transparent dots correspond to inclination angles of the system of 0 and 30°, respectively.

6. Conclusions

We have investigated the long-term spectral and temporal evolution of the ULX M 81 X-6, finding that the source oscillates between two different regimes, characterised by different values of hardness and luminosity. A modelling of the spectra in the two regimes indicates a probable high-energy component, which can be attributed to the emission coming from the accretion column, under the hypothesis of a NS accretor.

The two thermal components describing the broadband 0.3−10 keV spectra can be linked to the inner (higher temperature) and outer (lower temperature) regions of the accretion disc. The constant temperature of the softer of these components would suggest that the spectral changes are not due to variations in the mass-accretion rate. On the other hand, the amount of variation in luminosity is not high enough to justify transitions into or out of the propeller regime, assuming that M 81 X-6 has the same properties of other known ULXs. Therefore, the transitions between the two spectral states can be attributed to a precession of the accretion disc.

Assuming a NS accretor, we tested two possible scenarios: Lense-Thirring precession and precession due to the magnetic torque applied by the NS magnetic field onto the accretion disc. In the case of Lense-Thirring precession, we estimated the characteristic radii of the system, the magnetic field (B ≲ 10¹⁰ G), and the expected spin period (∼5 s) of the NS. In the case of a magnetic torque, a higher magnetic field (B ∼ 10¹² − 10¹⁴ G) is required, with a similar spin period. However, the latter mechanism would be disfavoured because it would require unrealistically high values of the magnetic field strength (B > 10¹⁵ G) for known PULXs.

Clearly, the direct detection of pulsations would be the strongest evidence in favour of the NS-ULX nature of M 81 X-6. However, a direct detection is difficult, not only because of instrumental limitations (e.g. time resolution, sensitivity, exposure times, etc.), but also because of the intrinsic physics and geometry of ULXs. For instance, in our case, the inferred low fraction of the flux of the cut-off power law component over the total flux, coupled with the upper limit on the PF of ∼10% found for one XMM-Newton/pn observation, would suggest a scenario where the pulsation is somehow ‘diluted’ by the wind cone in the outer photosphere. This may be a consequence of the small extent of the magnetosphere of the NS, compared to the size of the spherisation radius of the accretion disc, as inferred from the application of the Lense-Thirring precessing model (cf. Sect. 5.4).

New observational data and numerical simulations may prove crucial to more broadly understanding the properties of PULXs and thus discriminating between the two precession mechanisms. Moreover, future pointed longer observations of M 81 X-6 with XMM-Newton, NuSTAR, or next-generation X-ray satellites, such as XRISM (XRISM Science Team 2020) or Athena (Nandra et al. 2013), are encouraged to systematically look for pulsations and strengthen the findings of this work.

Acknowledgments

This research has made use of data obtained from the Chandra Data Archive and the Chandra Source Catalogue, and software provided by the Chandra X-ray Center (CXC) in the application package CIAO; of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Space Science Data Center (SSDC, Italy) and the California Institute of Technology (Caltech, USA); of observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This research has also made use of the software HENDRICS, based on the software package Stingray (Bachetti et al. 2022; Huppenkothen et al. 2019a,b). R.A., N.W., and O.G. acknowledge the support of the CNES for this work. The authors thank S. Vaughan for clarifications on the computation of the fractional variability.

References

All Tables

All Figures

Fig. 1. Temporal evolution of M 81 X-6. Left panel: Swift/XRT light curve in the 0.3−10 keV energy band. Each data point corresponds to an observation. The time interval spun by the Chandra observations considered in this work is coloured in light blue, while that of the most recent XMM-Newton observations is coloured in light orange or represented by a dotted orange vertical line. The vertical dotted green lines mark the oldest XMM-Newton observations (cf. G21) and the red one the NuSTAR observation used in this work. Right panel: histogram of Swift/XRT count rates, showing a clear bimodal distribution.

In the text

	Fig. 2. Hardness-intensity diagram (HID) of the Chandra (blue) and 2020 XMM-Newton (orange) datasets, in comparison with the old XMM-Newton observations (fainter text and points) as reported in G21.
In the text

	Fig. 3. Evolution of Swift/XRT broadband count rates (blue dots) as a function of the HR, in comparison with the latest Chandra (empty orange circles) and XMM-Newton data (empty orange diamonds). Data from different XMM-Newton-EPIC instruments are plotted individually. Errors on Chandra and XMM-Newton data are not shown.
In the text

	Fig. 4. Spectra of Chandra ObsID 18047 (blue) and XMM-Newton/pn ObsID 0870900201 (orange), representative of the HHL and SLL states, respectively. Spectra are fitted with the model tbabs*tbabs*(diskbb+diskbb) (solid line) with the two black-body components plotted as dashed lines. Bottom panel: residuals of the fit.
In the text

	Fig. 5. Light curves and HR of Chandra ObsID 18047. Top panel: light curves in the hard (orange) and soft (blue) energy bands, binned at 1000 s, with dashed lines representing the mean count rate of each light curve. Bottom panel: HR, computed as the number of counts in the hard band divided by the total number of counts. The dashed line corresponds to the mean HR.
In the text

Fig. 6. Bolometric luminosity as a function of the temperature for the soft (left panel) and hard (right panel) black-body components, in comparison with the older XMM-Newton observations as in G21 (in grey). For the soft component, three L − T laws are overplotted: L ∝ T4 for a thin disc (Shakura & Sunyaev 1973), L ∝ T2 for an advection-dominated disc around a BH (Watarai et al. 2000), and L ∝ T−4 as observed for NGC 1313 X-2 (Robba et al. 2021). The bolometric luminosities are computed in the 0.01−100 keV range.

In the text

	Fig. 7. HID of the Chandra observations and one XMM-Newton observation taken in 2001 (G21). The circle size is proportional to the fractional variability, only shown on observations for which it could be constrained.
In the text

Fig. 8. Lense-Thirring accretion plane adapted from Middleton et al. (2019). Left and right panels: resulting Tsph and precessing period of the wind for εwind = 0.2 and 0.8, respectively. The values expected for a 20 M⊙ BH with a dimensionless spin of 0.1 are shown in dashed dark blue and could be considered as a rough upper limit for the location of BHs in the plane. The rest of the lines show the values for NSs with different dipole magnetic field strengths, as indicated. The solid and dashed lines show the effect of varying the spin of the NS from 3 × 10−5 to 3 × 10−3, which we have converted to periods assuming canonical values for the moment of inertia and mass of the NS (see the text for details) and represent typical values observed in ULXs.

In the text

Fig. 9. Precessing period due to the interaction of the NS magnetic field with the accretion disc (Lipunov & Shakura 1980) as a function of the NS spin period calculated for different values of the NS dipole magnetic field. Left and right panels: calculations for ϵw = 0.2 and ϵw = 0.8, respectively. The horizontal dashed line indicates the super-orbital period of 110 days of M 81 X-6 (Lin et al. 2015), and we have considered the temperature of the soft component in the calculations of Rmag.

In the text

	Fig. 10. Comparison of the physical sizes of the system derived from the best-fit multi-temperature disc components (dots) and from the Lense-Thirring precession model (vertical lines; see Sect. 5.4 for details). The solid and transparent dots correspond to inclination angles of the system of 0 and 30°, respectively.
In the text