# The ultraluminous X-ray source M81 X-6: a weakly magnetised neutron star with a precessing accretion disc?

Roberta Amato, Andres Gúrpide, Natalie A. Webb, Olivier Godet, Matthew J. Middleton

<figure><img src="image_1.png" /></figure>

► To cite this version:

Roberta Amato, Andres Gúrpide, Natalie A. Webb, Olivier Godet, Matthew J. Middleton. The ultraluminous X-ray source M81 X-6: a weakly magnetised neutron star with a precessing accretion disc?. Astronomy and Astrophysics - A&A, 2023, 669, pp.A130. 10.1051/0004-6361/202244576. hal-03922288

HAL Id: hal-03922288

https://hal.science/hal-03922288

Submitted on 26 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# The ultraluminous X-ray source M81 X-6: a weakly magnetised neutron star with a precessing accretion disc?

R. Amato<fnref n="1" />, A. Gúrpide<fnref n="1" />,<fnref n="2" />, N. A. Webb<fnref n="1" />, O. Godet<fnref n="1" />, and M. J. Middleton<fnref n="2" />

Received 22 July 2022 / Accepted 10 December 2022

## ABSTRACT

**Context.** Thanks to their proximity, ultraluminous X-ray sources (ULXs) represent a privileged astrophysical laboratory to study super-Eddington accretion. Current open questions concern the nature of the compact object, which is still hard to determine in those cases where pulsations are not directly detected, and the mechanisms responsible for the spectral changes observed in many ULXs.

**Aims.** We investigate the nature of the ULX M 81 X-6, which has been suggested to harbour a neutron star (NS), by studying its long-term X-ray spectral and temporal evolution, with the goal of assessing the astrophysical phenomena responsible for its spectral changes.

**Methods.** Using the rich set of available archival data from *XMM-Newton*, *Chandra*, *NuSTAR*, and *Swift/XRT*, we tracked the evolution of the source on the hardness-intensity diagram and inferred the different emitting regions of the system and their geometry, as well as the mechanisms responsible for the spectral transitions.

**Results.** We find that the source oscillates between two main states: one characterised by a hard and luminous spectrum and the other at low hardness and luminosity. The properties of the soft component remain constant between the two states, suggesting that changes in the mass-transfer rate are not driving the spectral transitions. Instead, the bi-modal behaviour of the source and the known super-orbital period would point to the precession of the accretion disc. Here, we tested two theoretical models: (1) Lense-Thirring precession, which can explain the super-orbital period if the NS has a magnetic field $B \le 10^{10}$ G, supporting the idea of M 81 X-6 as a weakly magnetised NS, and (2) precession due to the torque of the NS magnetic field, which leads to $B \ge 10^{11}$ G. However, the latter scenario, assuming M 81 X-6 shares similar properties with other NS-ULXs, is disfavoured because it would require magnetic field strengths ($B > 10^{15}$ G) much higher than those known for other pulsating ULXs. We further show that the contribution from the hard component attributed to the putative accretion column sits just below the typical values found in pulsating ULXs, which, together with the low value of the pulsed fraction (≤10%) found for one *XMM-Newton/pn* observation, could explain the source's lack of pulsations.

**Conclusions.** The spectral properties and variability of M 81 X-6 can be accounted for if the accretor is a NS with a low magnetic field. Under the hypothesis of Lense-Thirring precession, we predict a spin period of the NS of a few seconds. We encourage future X-ray pointed observations to look for pulsations and/or spectral signatures of the magnetic field.

**Key words.** accretion, accretion disks – stars: neutron – X-rays: binaries – X-rays: individuals: M 81 X-6

## 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 \ge 10^{39}$ erg s$^{-1}$, 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^{41}$ erg s$^{-1}$), 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^{39}–10^{41}$ erg s$^{-1}$ 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

<footnote n="1">¹ IRAP, CNRS, Université de Toulouse, CNES, 9 Avenue du Colonel Roche, 31028 Toulouse, France
e-mail: roberta.amato@irap.omp.eu</footnote>

<footnote n="2">² School of Physics & Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>
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) \propto 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).

M81 X-6 is a bright ULX ($L_X \ge 2 \times 10^{39}$ erg s$^{-1}$) 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 $\ge 3 \times 10^{39}$ erg s$^{-1}$. 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 M81 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 M81 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 M81 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

M81 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 M81 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 M81 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$^{-1}$ 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).

<footnote n="1">¹ https://cxc.harvard.edu/ciao/ahelp/pileup_map.html</footnote>
<table><caption>**Table 1.** Log of *Chandra* and *XMM-Newton* observations considered in this work.</caption><thead><tr><th>Observatory</th><th>ObsID</th><th>Date</th><th>Instr.</th><th>Exp. (ks)</th><th>Net counts</th></tr></thead><tbody><tr><td rowspan="12">Chandra</td><td>18052</td><td>2016-02-01</td><td>ACIS-I</td><td>20.2</td><td>1492</td></tr><tr><td>18054</td><td>2016-06-21</td><td>ACIS-I</td><td>65.71</td><td>5825</td></tr><tr><td>18875</td><td>2016-06-24</td><td>ACIS-I</td><td>32.36</td><td>4236</td></tr><tr><td>18047</td><td>2016-07-04</td><td>ACIS-I</td><td>100.43</td><td>12570</td></tr><tr><td>18817</td><td>2016-08-07</td><td>ACIS-I</td><td>39.55</td><td>2978</td></tr><tr><td>19685</td><td>2016-08-08</td><td>ACIS-I</td><td>36.59</td><td>2264</td></tr><tr><td>18053</td><td>2017-01-05</td><td>ACIS-I</td><td>29.59</td><td>3296</td></tr><tr><td>19982</td><td>2017-01-12</td><td>ACIS-I</td><td>67.49</td><td>8510</td></tr><tr><td>18051</td><td>2017-01-24</td><td>ACIS-I</td><td>36.58</td><td>6504</td></tr><tr><td>19993</td><td>2017-01-26</td><td>ACIS-I</td><td>24.75</td><td>4116</td></tr><tr><td>19992</td><td>2017-01-28</td><td>ACIS-I</td><td>34.61</td><td>5967</td></tr><tr><td rowspan="10">XMM-Newton</td><td>0843840101</td><td>2020-03-18</td><td>MOS2</td><td>22.72</td><td>2078</td></tr><tr><td>0870930101</td><td>2020-10-17</td><td>MOS1/MOS2</td><td>22.27/22.28</td><td>1012/781</td></tr><tr><td>0870930301</td><td>2020-10-29</td><td>MOS2</td><td>15.65</td><td>503</td></tr><tr><td>0870930401</td><td>2020-11-06</td><td>MOS1/MOS2</td><td>5.43/4.9</td><td>159/118</td></tr><tr><td>0870900101</td><td>2020-11-14</td><td>pn/MOS2</td><td>17.15/24.45</td><td>4087/1815</td></tr><tr><td>0870900201</td><td>2020-11-16</td><td>pn/MOS1/MOS2</td><td>12.74/25.25/25.24</td><td>3854/2423/2544</td></tr><tr><td>0870930501</td><td>2020-11-18</td><td>MOS1</td><td>23.6</td><td>913</td></tr><tr><td>0870931001</td><td>2020-11-24</td><td>MOS1/MOS2</td><td>24.5/20.35</td><td>829/640</td></tr></tbody></table>

Notes. XMM-Newton exposure times have been cleaned by high-background intervals.

## 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<fnref n="2" />. For each of the 603 observations, we used the online tools<fnref n="3" /> (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.

## 2.4. NuSTAR

*NuSTAR* has observed the galaxy M81 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 M81. 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, M81 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 \times 10^{20} \text{ cm}^{-2}$ (HI4PI Collaboration 2016) and local absorption was fixed to $11.7 \times 10^{20} \text{ cm}^{-2}$, 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 M81 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 $\chi^2$ 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 $\sim 2–5 \times 10^{39} \text{ erg s}^{-1}$.

We noticed that in almost all *Chandra* observations and in two *XMM-Newton* ones (ObsID 0870900101 and 0870900201) residuals at energies $\le 1$ keV are marginally evident. Those residuals have sometimes been interpreted as signatures of powerful outflows irradiated from the surface of the accretion disc

<footnote n="2">² https://www.swift.ac.uk/swift.live/index.php</footnote>

<footnote n="3">³ https://www.swift.ac.uk/user_objects/</footnote>
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 2.** Best-fit parameters for each *Chandra* observation.</caption><thead><tr><th>ObsID</th><th>kT<sub>soft</sub> (keV)</th><th>Norm.</th><th>kT<sub>hard</sub> (keV)</th><th>Norm. (10<sup>-2</sup>)</th><th>Flux (10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Lum. (10<sup>39</sup> erg s<sup>-1</sup>)</th><th>χ<sup>2</sup><sub>red</sub>/d.o.f.</th><th>State</th></tr></thead><tbody><tr><td>G21</td><td>0.4–0.6</td><td>0.3–2</td><td>1.4–2.5</td><td>0.1–3</td><td></td><td>2.4–7.4</td><td>0.6–1.2</td><td></td></tr><tr><td>18052</td><td>0.60 ± 0.15</td><td>0.28<sup>+0.43</sup><sub>-0.14</sub></td><td>1.7<sup>(b)</sup></td><td>0.07<sup>+1.3</sup><sub>-3.4</sub></td><td>1.44 ± 0.06</td><td>2.2 ± 0.1</td><td>0.72/78</td><td>SLL</td></tr><tr><td>18054</td><td>0.73 ± 0.09</td><td>0.19<sup>+0.08</sup><sub>-0.05</sub></td><td>3<sup>(b)</sup></td><td>0.02<sup>+0.09</sup><sub>-0.02</sub></td><td>2.14 ± 0.05</td><td>3.34 ± 0.07</td><td>0.90/206</td><td>SLL</td></tr><tr><td>18875</td><td>0.8 ± 0.3</td><td>0.14<sup>+0.36</sup><sub>-0.07</sub></td><td>1.6<sup>(b)</sup></td><td>0.4<sup>+1.2</sup><sub>-0.4</sub></td><td>2.89 ± 0.07</td><td>4.15 ± 0.11</td><td>0.93/182</td><td>HHL</td></tr><tr><td>18047</td><td>0.35<sup>+0.12</sup><sub>-0.08</sub></td><td>2.28<sup>+5.6</sup><sub>-1.6</sub></td><td>1.65<sup>+0.14</sup><sub>-0.09</sub></td><td>1.8 ± 0.5</td><td>3.26 ± 0.05</td><td>5.08 ± 0.08</td><td>0.98/321</td><td>HHL</td></tr><tr><td>18817</td><td>0.46<sup>+0.15</sup><sub>-0.12</sub></td><td>0.8<sup>+1.7</sup><sub>-0.5</sub></td><td>2.0<sup>+1.8</sup><sub>-0.5</sub></td><td>0.3<sup>+0.6</sup><sub>-0.3</sub></td><td>1.69 ± 0.05</td><td>2.64 ± 0.08</td><td>0.93/140</td><td>SLL</td></tr><tr><td>19685</td><td>0.40<sup>+0.09</sup><sub>-0.07</sub></td><td>1.5<sup>+2.1</sup><sub>-0.9</sub></td><td>1.8<sup>+0.5</sup><sub>-0.3</sub></td><td>0.5<sup>+0.5</sup><sub>-0.3</sub></td><td>1.66 ± 0.05</td><td>2.59 ± 0.09</td><td>1.01/129</td><td>SLL</td></tr><tr><td>18053</td><td>0.4<sup>+5.1</sup><sub>-0.3</sub></td><td>0.7<sup>+12</sup><sub>-0.6</sub></td><td>1.3<sup>(b)</sup></td><td>3.4<sup>+1.2</sup><sub>-3.4</sub></td><td>2.80 ± 0.08</td><td>4.36 ± 0.13</td><td>0.98/154</td><td>HHL</td></tr><tr><td>19982</td><td>0.32<sup>+0.56</sup><sub>-0.17</sub></td><td>1.18<sup>+20</sup><sub>-1.11</sub></td><td>1.45<sup>+0.27</sup><sub>-0.07</sub></td><td>3.3<sup>+0.7</sup><sub>-2.0</sub></td><td>3.17 ± 0.06</td><td>4.94 ± 0.09</td><td>0.97/274</td><td>HHL</td></tr><tr><td>18051</td><td>0.6<sup>+3.3</sup><sub>-0.6</sub></td><td>0.1<sup>(a)</sup></td><td>1.4<sup>+4.0</sup><sub>-0.1</sub></td><td>3.4<sup>+1.3</sup><sub>-3.4</sub></td><td>3.25 ± 0.07</td><td>5.1 ± 0.1</td><td>1.13/232</td><td>HHL</td></tr><tr><td>19993</td><td>1.17<sup>(a)</sup></td><td>0.06<sup>(a)</sup></td><td>2.8<sup>(a)</sup></td><td>0.075<sup>(a)</sup></td><td>3.08 ± 0.08</td><td>4.81 ± 0.13</td><td>0.84/180</td><td>HHL</td></tr><tr><td>19992</td><td>0.74<sup>(a)</sup></td><td>0.06<sup>(a)</sup></td><td>1.4<sup>(b)</sup></td><td>2.4<sup>(a)</sup></td><td>3.04 ± 0.07</td><td>4.74 ± 0.11</td><td>1.01/216</td><td>HHL</td></tr></tbody></table>

**Notes.** The interstellar absorption is fixed to $9.8 \times 10^{20} \text{ cm}^{-2}$ and the local one to $11.7 \times 10^{20} \text{ cm}^{-2}$ (G21). Errors are at 90%. Fluxes and luminosities are both unabsorbed and computed in the 0.3–10 keV range. The range of the best-fit values from G21 are reported in the first row for an easy reference to the reader. The last column reports whether the source was caught in the hard/high-luminosity (HHL) or soft/low-luminosity (SLL) regime. $^{(a)}$Parameter not constrained. $^{(b)}$Lower limit.

<table><caption>**Table 3.** Best-fit parameters for each *XMM-Newton* observation.</caption><thead><tr><th>ObsID</th><th>kT<sub>soft</sub> (keV)</th><th>Norm.</th><th>kT<sub>hard</sub> (keV)</th><th>Norm. (10<sup>-2</sup>)</th><th>Flux (10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Lum. (10<sup>39</sup> erg s<sup>-1</sup>)</th><th>χ<sup>2</sup><sub>red</sub>/d.o.f.</th><th>State</th></tr></thead><tbody><tr><td>G21</td><td>0.4–0.6</td><td>0.3–2</td><td>1.4–2.5</td><td>0.1–3</td><td></td><td>2.4–7.4</td><td>0.6–1.2</td><td></td></tr><tr><td>0843840101</td><td>0.43<sup>+0.16</sup><sub>-0.12</sub></td><td>0.8<sup>+1.6</sup><sub>-0.5</sub></td><td>1.7<sup>+0.4</sup><sub>-0.3</sub></td><td>1.4<sup>+1.2</sup><sub>-0.8</sub></td><td>2.95 ± 0.11</td><td>4.59 ± 0.17</td><td>0.97/85</td><td>HHL</td></tr><tr><td>0870930101</td><td>0.43<sup>+0.08</sup><sub>-0.07</sub></td><td>1.1<sup>+1.0</sup><sub>-0.5</sub></td><td>2.4<sup>+2.4</sup><sub>-0.7</sub></td><td>0.17<sup>+0.5</sup><sub>-0.16</sub></td><td>1.79 ± 0.07</td><td>2.79 ± 0.11</td><td>0.86/79</td><td>SLL</td></tr><tr><td>0870930301</td><td>0.33<sup>+0.12</sup><sub>-0.09</sub></td><td>2.2<sup>+5.3</sup><sub>-1.5</sub></td><td>2.3<sup>+3.5</sup><sub>-0.8</sub></td><td>0.2<sup>+0.7</sup><sub>-0.2</sub></td><td>1.58 ± 0.13</td><td>2.5 ± 0.2</td><td>1.11/23</td><td>SLL</td></tr><tr><td>0870930401</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>SSL</td></tr><tr><td>0870900101</td><td>0.45<sup>+0.06</sup><sub>-0.05</sub></td><td>0.9<sup>+0.5</sup><sub>-0.3</sub></td><td>1.9<sup>+0.4</sup><sub>-0.2</sub></td><td>0.6<sup>+0.5</sup><sub>-0.3</sub></td><td>2.30 ± 0.05</td><td>3.59 ± 0.08</td><td>1.15/150</td><td>SLL</td></tr><tr><td>0870900201</td><td>0.42 ± 0.03</td><td>1.2<sup>+0.4</sup><sub>-0.3</sub></td><td>2.0<sup>+0.3</sup><sub>-0.2</sub></td><td>0.32<sup>+0.21</sup><sub>-0.14</sub></td><td>1.64 ± 0.03</td><td>2.57 ± 0.05</td><td>1.15/233</td><td>SLL</td></tr><tr><td>0870930501</td><td>0.43 ± 0.11</td><td>1.0<sup>+1.7</sup><sub>-0.6</sub></td><td>2.4<sup>+6.9</sup><sub>-0.9</sub></td><td>0.17<sup>+0.9</sup><sub>-0.17</sub></td><td>1.6 ± 0.1</td><td>2.57 ± 0.15</td><td>0.77/45</td><td>SLL</td></tr><tr><td>0870931001</td><td>0.40<sup>+0.06</sup><sub>-0.05</sub></td><td>1.4<sup>+1.1</sup><sub>-0.6</sub></td><td>4.5<sup>+4.5</sup><sub>-2.1</sub></td><td>0.02<sup>+0.14</sup><sub>-0.02</sub></td><td>1.75 ± 0.08</td><td>2.73 ± 0.13</td><td>0.84/78</td><td>SLL</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

(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<sup>39</sup> erg s<sup>-1</sup>), while in the recent observations it never went above L<sub>X</sub> ~ 5 × 10<sup>39</sup> erg s<sup>-1</sup>.

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).

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<sup>39</sup> erg s<sup>-1</sup>, and a soft, low-luminosity (SLL) state, with HR ≲ 2 and L<sub>X</sub> ~ 3 × 10<sup>39</sup> erg s<sup>-1</sup> (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).

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 T<sub>hard</sub> = 1.46<sup>+0.17</sup><sub>-0.05</sub> keV and T<sub>soft</sub> = 0.45<sup>+0.51</sup><sub>-0.15</sub> keV for the HHL regime, and T<sub>hard</sub> = 2.2<sup>+0.3</sup><sub>-0.2</sub> keV and T<sub>soft</sub> = 0.48<sup>+0.05</sup><sub>-0.04</sub> 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 n<sub>H</sub> = (11.2<sup>+8.1</sup><sub>-5.0</sub>) × 10<sup>20</sup> cm<sup>-2</sup> and n<sub>H</sub> = (7.7<sup>+2.3</sup><sub>-2.1</sub>) × 10<sup>20</sup> cm<sup>-2</sup>, for the HHL and the SLL regimes, respectively, consistent in both cases with the average value of 11.7 × 10<sup>20</sup> cm<sup>-2</sup> 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
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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 M81 X-6. To constrain the parameters of the cut-off power law, we used the *NuSTAR* observation, taken when M81 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 Ecut 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
<table><caption>**Table 4.** Best-fit parameters of the observations in the two regimes.</caption><thead><tr><th></th><th>HHL</th><th>SLL</th><th>HHL<br>+cutoffpl</th><th>SLL<br>+cutoffpl</th></tr></thead><tbody><tr><td>N<sub>H</sub> (10<sup>20</sup> cm<sup>-2</sup>)</td><td>11.2<sup>+8.1</sup><sub>-5.0</sub></td><td>7.7<sup>+2.3</sup><sub>-2.1</sub></td><td>11.5<sup>+7.9</sup><sub>-5.0</sub></td><td>9.40 ± 0.03</td></tr><tr><td>kT<sub>soft</sub> (keV)</td><td>0.45<sup>+0.51</sup><sub>-0.15</sub></td><td>0.48<sup>+0.05</sup><sub>-0.04</sub></td><td>0.4<sup>+0.5</sup><sub>-0.2</sub></td><td>0.39 ± 0.09</td></tr><tr><td>Norm.</td><td>0.3<sup>+2.3</sup><sub>-0.2</sub></td><td>0.6<sup>+0.3</sup><sub>-0.2</sub></td><td>0.3<sup>+2.6</sup><sub>-0.3</sub></td><td>1.2<sup>+1.9</sup><sub>-0.5</sub></td></tr><tr><td>kT<sub>hard</sub> (keV)</td><td>1.46<sup>+0.17</sup><sub>-0.05</sub></td><td>2.2<sup>+0.3</sup><sub>-0.2</sub></td><td>1.46<sup>+0.16</sup><sub>-0.15</sub></td><td>0.9<sup>+0.8</sup><sub>-0.2</sub></td></tr><tr><td>Norm. (10<sup>-2</sup>)</td><td>3.1<sup>+0.5</sup><sub>-1.5</sub></td><td>0.22<sup>+0.13</sup><sub>-0.09</sub></td><td>3.0<sup>+1.4</sup><sub>-0.9</sub></td><td>3.6<sup>+13.6</sup><sub>-3.3</sub></td></tr><tr><td>Γ</td><td>—</td><td>—</td><td>0.38</td><td>0.38</td></tr><tr><td>E<sub>cut</sub> (keV)</td><td>—</td><td>—</td><td>4.3</td><td>4.3</td></tr><tr><td>Norm. (10<sup>-5</sup>)</td><td>—</td><td>—</td><td>2.8 <sup>(a)</sup></td><td>6.8<sup>+0.7</sup><sub>-2.6</sub></td></tr><tr><td>F<sub>tot</sub> (10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</td><td>3.06 ± 0.03</td><td>1.57 ± 0.02</td><td>3.06 ± 0.03</td><td>1.64 ± 0.02</td></tr><tr><td>F<sub>cutoff</sub>/F<sub>tot</sub></td><td>—</td><td>—</td><td>0.013 <sup>(a)</sup></td><td>0.50 ± 0.02</td></tr><tr><td>L<sub>X</sub> (10<sup>39</sup> erg s<sup>-1</sup>)</td><td>4.77 ± 0.05</td><td>2.44 ± 0.03</td><td>4.77 ± 0.05</td><td>2.56 ± 0.03</td></tr><tr><td>χ<sup>2</sup>/d.o.f.</td><td>1220.12/1153</td><td>884.20/832</td><td>1220.11/1152</td><td>880.08/831</td></tr><tr><td>n.h.p. <sup>(b)</sup></td><td>0.083</td><td>0.102</td><td>0.08</td><td>0.116</td></tr></tbody></table>

**Notes.** The unabsorbed X-ray fluxes and luminosities are in the 0.3–10 keV range, while the flux ratio $F_{cutoff}/F_{tot}$ is computed in the range 0.3–40 keV, for the purpose of comparison with Walton et al. (2018a). The cut-off power law parameters, except the normalisation, are frozen to the best-fit ones of the SLL regime obtained using the NuSTAR dataset; the interstellar absorption is fixed to $N_H = 9.8 \times 10^{20} \text{ cm}^{-2}$ (G21). Errors are at 90% confidence level. <sup>(a)</sup>Upper limit. <sup>(b)</sup>The null hypothesis probability (n.h.p.) is the probability of the observed data being drawn from the model, for a certain value of the $\chi^2$ and the degrees of freedom (d.o.f.). In the present work, we chose the common threshold of 0.05 and, hence, a fit with n.h.p. >0.05 can be considered statistically acceptable.

be $E_{\text{cut}} = 4.3_{-0.3}^{+1.3}$ keV and for the power law index $\Gamma = 0.38_{-0.90}^{+0.15}$. 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_{\text{soft}} \approx 0.4$ keV), while the best-fit temperature of the hard black body seems to be higher in the HHL state ($T_{\text{hard}} = 1.46_{-0.15}^{+0.16}$ keV) and lower in the SLL state ($T_{\text{hard}} = 0.9_{-0.2}^{+0.8}$ 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 $\chi^2_{\text{red}}(\text{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 $T_{\text{soft}} = 0.34_{-0.05}^{+0.06}$ keV, with a normalisation of $1.4_{-0.8}^{+1.7}$. 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_{\text{soft}} \sim 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 \propto T^4$ in the case of a thin disc (Shakura & Sunyaev 1973) and $L \propto T^2$ 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 \propto T^{-4}$ 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.

As mentioned in Sect. 1, the HHL state of M 81 X-6 has sometimes been also modelled with a BD component
<figure><img src="image_8.png" /><figcaption>**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* ~ *T*⁴ for a thin disc (Shakura & Sunyaev 1973), *L* ~ *T*² for an advection-dominated disc around a BH (Watarai et al. 2000), and *L* ~ *T*⁻⁴ as observed for NGC 1313 X-2 (Robba et al. 2021). The bolometric luminosities are computed in the 0.01–100 keV range.</figcaption></figure>

(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 *kT* = 1.50<sup>+0.07</sup><sub>-0.19</sub> keV, a parameter *p* = 0.70<sup>+0.07</sup><sub>-0.03</sub>, and a normalisation of 2.4<sup>+2.0</sup><sub>-0.6</sub> × 10<sup>-2</sup>, with a *χ²*<sub>red</sub>(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

<table><caption>**Table 5.** Fractional variability of some of the *Chandra* observations. *XMM-Newton* ObsID 01118001001 was taken in 2001, and it is analysed in G21.</caption><thead><tr><td>ObsID</td><td>F<sub>var</sub></td><td>N. segments</td><td>State</td></tr></thead><tbody><tr><td>18047</td><td>0.09 ± 0.02</td><td>3</td><td>HHL</td></tr><tr><td>18817</td><td>0.04 ± 0.06</td><td>2</td><td>SLL</td></tr><tr><td>18875</td><td>0.07 ± 0.03</td><td>1</td><td>HHL</td></tr><tr><td>19685</td><td>0.04 ± 0.07</td><td>1</td><td>SLL</td></tr><tr><td>19982</td><td>0.04 ± 0.04</td><td>1</td><td>HHL</td></tr><tr><td>0111800101</td><td>0.037 ± 0.005</td><td>4</td><td>HHL</td></tr></tbody></table>

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).

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 $Z_n^2$ 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*–*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
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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. M81 X-6 spectral states**

According to the general classification of ULXs (Sutton et al. 2013), with its HR always $\ge 1$ (see Figs. 2 and 3), M81 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 $\sim$2–3 and $L_X \sim 4–6 \times 10^{39}$ erg s$^{-1}$) and SLL states (HR $\le 2$ and $L_X \sim 3 \times 10^{39}$ erg s$^{-1}$; 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 $\sim 1–3$) and luminosity ($L_X \sim 2–10 \times 10^{39}$ erg s$^{-1}$) 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 $\ge 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. M81 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 $\le 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} \le 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 ≈170 $P^{2/3} M_{1.4}^{1/3} R_6^{-1}$ (Tsygankov et al. 2016, Eq. (5))<fnref n="4" />, where P is the NS rotational period, $M_{1.4}$ its mass in units of $1.4 M_\odot$, and $R_6$ its radius in units of 10 km. Using the canonical values for a NS of $1.4 M_\odot$ and 10-km radius, a quick estimation would lead to $P \sim 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_{\text{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_{\text{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_{\text{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, $\beta$) 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, $\varepsilon_{\text{wind}} = L_{\text{wind}}/(L_{\text{wind}} + L_{\text{rad}})$, with $L_{\text{wind}}$ the kinetic luminosity of the wind and $L_{\text{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^{45}$ g cm²) and its mass ($M = 1.4 M_\odot$), we find that $3 \times 10^{-2} > a_* > 3 \times 10^{-5}$ should approximately cover the range of observed periods in ULXs (~0.1 s < $P$ < 10 s). For $\varepsilon_{\text{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_\odot$ 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_{\text{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_{\text{sph}}$. Because of the latter reason, our naive expectation would be that our estimate for $T_{\text{sph}}$ constitutes an upper limit. These caveats should be borne in

<footnote n="4">⁴ Note that this formulation assumes sub-Eddington accretion and might be erroneous for ULXs. A correction of this formula is beyond the goals of this paper (see Middleton et al. 2023).</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Lense-Thirring accretion plane adapted from Middleton et al. (2019). Left and right panels: resulting $T_{\text{sph}}$ and precessing period of the wind for $\epsilon_{\text{wind}} = 0.2$ and 0.8, respectively. The values expected for a $20 M_{\odot}$ 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 \times 10^{-5}$ to $3 \times 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.</figcaption></figure>

mind when interpreting these results and therefore our estimates should be considered semi-quantitative. Improvement in theoretical modelling of PULXs should help reduce these uncertainties in the future.

From Fig. 8, we see that the position of M 81 X-6 on the plane indicates a weakly magnetised NS ($B < 10^{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^{10}$ G. This is because for low magnetic field values ($B < 10^{10}$ G) and accretion rates with $T_{\text{sph}} \lesssim 0.6$ keV we find already that $R_{\text{mag}} \rightarrow R_{\text{ISCO}}$ (i.e. the magnetospheric radius $R_{\text{mag}}$ is pushed all the way to the innermost stable orbit $R_{\text{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^{12}$ G if we consider the curves for $\epsilon_w = 0.2$. We note however that mostly $B < 10^{10}$ G would seem favoured as $\epsilon_w = 0.2$ constitutes a limiting case and high values would rule out the case for $B \sim 10^{12}$ G. Moreover, since higher values of $\epsilon_w$ have the effect of increasing the size of the magnetosphere, as more mass is lost into the outflow, for $\epsilon_w > 0.3$, a magnetic field with $B = 10^{12}$ 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_{\text{sph}} \sim 0.4$ keV, $\epsilon_{\text{wind}} = 0.2$ and using $B = 10^{10}$ G, we would expect for M 81 X-6 a spin period $P_{\text{NS}} \sim 4.5$ s, a mass-accretion rate at the magnetospheric radius $\dot{m}_0 = 42$, a mass-accretion rate at the inner disc radius $\dot{m}_{\text{in}} \sim 37$, a magnetospheric radius $R_{\text{mag}} \sim 13 R_{\text{g}}$, a spherisation radius $R_{\text{sph}} \sim 302 R_{\text{g}}$, and an outer atmosphere radius $R_{\text{out}} \sim 758 R_{\text{g}}$, where $R_{\text{g}}$ is the gravitational radius ($R_{\text{g}} = GM_{\text{NS}}/c^2$). For $\epsilon_{\text{wind}} = 0.8$, we found $P_{\text{NS}} \sim 4.7$ s, $\dot{m}_0 = 21$, $\dot{m}_{\text{in}} \sim 12$, $R_{\text{mag}} \sim 18 R_{\text{g}}$, $R_{\text{sph}} \sim 129 R_{\text{g}}$, and $R_{\text{out}} \sim 726 R_{\text{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 $\dot{m}_0$ and $B$ found for both $\epsilon_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 rea-

sonable 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_{\text{sph}} > R_{\text{mag}}$, as suggested by G21. The overall low $R_{\text{mag}}$, close to the $R_{\text{ISCO}} = 6 R_{\text{g}}$, coupled with the large $R_{\text{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_{\text{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 $\ge 10^{11}$ 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)

$$P_{\text{mag,prec}}^{\text{mag}} = 1.5 \times 10^4 \left(\frac{B}{10^{12} \text{ G}}\right)^{-2} \left(\frac{R_{\text{NS}}}{10^6 \text{ cm}}\right)^{-2} \left(\frac{R_{\text{mag}}}{10^8 \text{ cm}}\right)^3 \\ \times \left(\frac{P_{\text{NS}}}{1 \text{ s}}\right)^{-1} \left(\frac{I}{10^{45} \text{ gr cm}^2}\right) \frac{1}{\cos \psi(3 \cos^2 \xi - 1)} \text{yr,} \quad (1)$$

where $\psi$ is the angle between the NS rotational axis and the normal to the accretion disc and $\xi$ is the angle between the NS magnetic field and its rotational axis. $R_{\text{NS}}$ is the radius of the NS at its surface, assumed to be 10 km. As above, we solve numerically for $R_{\text{mag}}$, based on $\dot{m}_0$ estimated from the temperature of the soft component and for different NS magnetic fields. We set $\xi = \psi = 30^{\circ}$ for reference and carry out the calculations for $\epsilon_w = 0.2$ and 0.8 as above and for $P_{\text{NS}} = 0.1-20$ s. We present the calculations in Fig. 9, which can also be useful for future reference.
<figure><img src="image_11.png" /></figure>

**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 $\epsilon_w = 0.2$ and $\epsilon_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 $R_{\text{mag}}$.

<figure><img src="image_12.png" /></figure>

**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.

As opposed to the Lense-Thirring scenario, this model would require a strongly magnetised ($10^{12}$ G $< B < 10^{14}$ 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^{15}$ 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 $\le 10^{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_{\text{sph}}$ and $R_{\text{out}}$, whereas the harder emission would come from the inner region from within $R_{\text{sph}}$ down to $R_{\text{ISCO}}$ for a BH-ULX or to $R_{\text{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_{\text{sph}}$ or fall within $R_{\text{sph}}$ and $R_{\text{out}}$. The only exception may be the case in the HHL state with $\epsilon_{\text{wind}} = 0.2$ (Fig. 10, top left panel), where the apparent inner radius falls slightly before $R_{\text{sph}}$, but, in any case, the soft-emission region would still develop in between the $R_{\text{sph}}$ and the $R_{\text{out}}$. On the other hand, we find relatively good agreement between $R_{\text{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_{\text{mag}}$ and $R_{\text{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.

## 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 temper-
ature) 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 vari-
ation 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 M81 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 preces-
sion 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.

Acknowledgements. 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

Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, ASP Conf. Ser., 101, 17

Bachetti, M. 2018, Astrophysics Source Code Library [record ascl:1805.019]

Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202

Bachetti, M., Huppenkothen, D., Khan, U., et al. 2022, https://doi.org/10.5281/zenodo.6394742

Brightman, M., Bachetti, M., Earnshaw, H., et al. 2022, ApJ, 925, 18

Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165

Carpano, S., Haberl, F., Maitra, C., & Vasilopoulos, G. 2018, MNRAS, 476, L45

Cseh, D., Corbel, S., Kaaret, P., et al. 2012, ApJ, 749, 17

Earnshaw, H. P., Roberts, T. P., & Sathyaprakash, R. 2018, MNRAS, 476, 4272

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177

Fabbiano, G. 1988, ApJ, 325, 544

Farrell, S. A., Webb, N. A., Barret, D., Godet, O., & Rodrigues, J. M. 2009, Nature, 460, 73

Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668, 417

Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, eds. D. R. Silva, & R. E. Doxsey, 6270, 62701V

Fürst, F., Walton, D. J., Harrison, F. A., et al. 2016, ApJ, 831, L14

Gao, Y., Wang, Q. D., Appleton, P. N., & Lucas, R. A. 2003, ApJ, 596, L171

Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005

Gladstone, J. C., Roberts, T. P., & Done, C. 2009, MNRAS, 397, 1836

Gúrpide, A., Godet, O., Koliopanos, F., Webb, N., & Olive, J. F. 2021a, A&A, 649, A104

Gúrpide, A., Godet, O., Vasilopoulos, G., Webb, N. A., & Olive, J. F. 2021b, A&A, 654, A10

Gúrpide, A., Parra, M., Godet, O., Contini, T., & Olive, J. F. 2022, A&A, 666, A100

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

HI4PI Collaboration (Ben Bekhti, N., et al.) 2016, A&A, 594, A116

Huppenkothen, D., Bachetti, M., Stevens, A., et al. 2019a, J. Open Sour. Softw., 4, 1393

Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019b, ApJ, 881, 39

Israel, G. L., Belfiore, A., Stella, L., et al. 2017a, Science, 355, 817

Israel, G. L., Papitto, A., Esposito, P., et al. 2017b, MNRAS, 466, L48

Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1

Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106

Jithesh, V., & Misra, R. 2018, ApJ, 868, 125

Jithesh, V., Anjana, C., & Misra, R. 2020, MNRAS, 494, 4026

Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303

King, A., & Lasota, J.-P. 2016, MNRAS, 458, L10

Koliopanos, F., Vasilopoulos, G., Godet, O., et al. 2017, A&A, 608, A47

Kosec, P., Pinto, C., Fabian, A. C., & Walton, D. J. 2018, MNRAS, 473, 5680

Kosec, P., Pinto, C., Reynolds, C. S., et al. 2021, MNRAS, 508, 3569

Kubota, A., Tanaka, Y., Makishima, K., et al. 1998, PASJ, 50, 667

Lin, L. C.-C., Hu, C.-P., Kong, A. K. H., et al. 2015, MNRAS, 454, 1644

Lipunov, V. M., & Shakura, N. I. 1980, Soviet Astron. Lett., 6, 14

Makishima, K., Maejima, Y., Mitsuda, K., et al. 1986, ApJ, 308, 635

Mezcua, M. 2017, Int. J. Mod. Phys. D, 26, 1730021

Middleton, M. J., & King, A. 2017, MNRAS, 470, L69

Middleton, M. J., Roberts, T. P., Done, C., & Jackson, F. E. 2011, MNRAS, 411, 644

Middleton, M. J., Heil, L., Pintore, F., Walton, D. J., & Roberts, T. P. 2015a, MNRAS, 447, 3243

Middleton, M. J., Walton, D. J., Fabian, A., et al. 2015b, MNRAS, 454, 3134

Middleton, M. J., Walton, D. J., Fabian, A., et al. 2018, MNRAS, 475, 154

Middleton, M. J., Fragile, P. C., Ingram, A., & Roberts, T. P. 2019, MNRAS, 489, 282

Middleton, M., Gúrpide, A., & Walton, D. J. 2023, MNRAS, 519, 2224

Mineshige, S., Hirano, A., Kitamoto, S., Yamada, T. T., & Fukue, J. 1994, ApJ, 426, 308

Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741

Mitsuda, K., Bautz, M., Inoue, H., et al. 2007, PASJ, 59, S1

Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Ingram, A. 2017, MNRAS, 467, 1202

Mushtukov, A. A., Ingram, A., Middleton, M., Nagirner, D. I., & van der Klis, M. 2019, MNRAS, 484, 687

Mushtukov, A. A., Portegies Zwart, S., Tsygankov, S. S., Nagirner, D. I., & Poutanen, J. 2021, MNRAS, 501, 2424

Nandra, K., Barret, D., Barcons, X., et al. 2013, ArXiv e-prints
[arXiv:1306.2307]

Pakull, M. W., & Mirioni, L. 2002, ApJ, submitted [arXiv:astro-ph/0202488]

Pinto, C., Middleton, M. J., & Fabian, A. C. 2016, Nature, 533, 64

Pinto, C., Walton, D. J., Kara, E., et al. 2020, MNRAS, 492, 4646

Pinto, C., Soria, R., Walton, D. J., et al. 2021, MNRAS, 505, 5058

Pintore, F., Zampieri, L., Stella, L., et al. 2017, ApJ, 836, 113
Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187

Quintin, E., Webb, N. A., Gúrpide, A., Bachetti, M., & Fürst, F. 2021, MNRAS, 503, 5485

Robba, A., Pinto, C., Walton, D. J., et al. 2021, A&A, 652, A118

Rodríguez Castillo, G. A., Israel, G. L., Belfiore, A., et al. 2020, ApJ, 895, 60

Salvaggio, C., Wolter, A., Pintore, F., et al. 2022, MNRAS, 512, 1814

Sathyaprakash, R., Roberts, T. P., Walton, D. J., et al. 2019, MNRAS, 488, L35

Sathyaprakash, R., Roberts, T. P., Grisé, F., et al. 2022, MNRAS, 511, 5346

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Song, X., Walton, D. J., Lansbury, G. B., et al. 2020, MNRAS, 491, 1260

Soria, R., Pakull, M. W., Motch, C., et al. 2021, MNRAS, 501, 1644

Sutton, A. D., Roberts, T. P., & Middleton, M. J. 2013, MNRAS, 435, 1758

Swartz, D. A., Ghosh, K. K., McCollough, M. L., et al. 2003, ApJS, 144, 213

Takeuchi, S., Ohsuga, K., & Mineshige, S. 2013, PASJ, 65, 88

Tsygankov, S. S., Mushtukov, A. A., Suleimanov, V. F., & Poutanen, J. 2016, MNRAS, 457, 1101

Urquhart, R., & Soria, R. 2016, MNRAS, 456, 1859

Vasilopoulos, G., Petropoulou, M., Koliopanos, F., et al. 2019, MNRAS, 488, 5225

Vasilopoulos, G., Lander, S. K., Koliopanos, F., & Bailyn, C. D. 2020, MNRAS, 491, 4949

Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Walton, D. J., Middleton, M. J., Pinto, C., et al. 2016a, ApJ, 826, L26

Walton, D. J., Fürst, F., Bachetti, M., et al. 2016b, ApJ, 827, L13

Walton, D. J., Fürst, F., Heida, M., et al. 2018a, ApJ, 856, 128

Walton, D. J., Fürst, F., Harrison, F. A., et al. 2018b, MNRAS, 473, 4360

Walton, D. J., Pinto, C., Nowak, M., et al. 2020, MNRAS, 494, 6012

Watarai, K.-Y., Fukue, J., Takeuchi, M., & Mineshige, S. 2000, PASJ, 52, 133

Weisskopf, M. C., Tananbaum, H. D., Van Speybroeck, L. P., & O'Dell, S. L. 2000, in X-Ray Optics, Instruments, and Missions III, eds. J. E. Truemper, & B. Aschenbach, SPIE Conf. Ser., 4012, 2

Weng, S.-S., & Feng, H. 2018, ApJ, 853, 115

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

XRISM Science Team 2020, ApJ, submitted [arXiv:2003.04962]