# Orbital Decay in M82 X-2

Matteo Bachetti, Marianne Heida, Thomas Maccarone, Daniela Huppenkothen, Gian Luca Israel, Didier Barret, Murray Brightman, Mckinley Brumback, Hannah P. Earnshaw, Karl Forster, et al.

► To cite this version:

Matteo Bachetti, Marianne Heida, Thomas Maccarone, Daniela Huppenkothen, Gian Luca Israel, et al. Orbital Decay in M82 X-2. The Astrophysical Journal, 2022, 937, 10.3847/1538-4357/ac8d67. insu-03867383

HAL Id: insu-03867383

https://insu.hal.science/insu-03867383

Submitted on 23 Nov 2022

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_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

Orbital Decay in M82 X-2

Matteo Bachetti¹ ID, Marianne Heida² ID, Thomas Maccarone³ ID, Daniela Huppenkothen⁴ ID, Gian Luca Israel⁵ ID, Didier Barret⁶ ID,
Murray Brightman⁷ ID, McKinley Brumback⁷ ID, Hannah P. Earnshaw⁷ ID, Karl Forster⁷ ID, Felix Fürst⁸ ID,
Brian W. Grefenstette⁷ ID, Fiona A. Harrison⁷, Amruta D. Jaodand⁷ ID, Kristin K. Madsen⁹ ID, Matthew Middleton¹⁰ ID,
Sean N. Pike⁷ ID, Maura Pilia¹ ID, Juri Poutanen¹¹,¹² ID, Daniel Stern⁷ ID, John A. Tomsick¹³ ID, Dominic J. Walton¹⁴,¹⁵ ID,
Natalie Webb⁶, and Jörn Wilms¹⁶ ID

¹ INAF-Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy; matteo.bachetti@inaf.it

² European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany

³Department of Physics and Astronomy, Texas Tech University, Lubbock, TX 79409-1051, USA

⁴SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA, Utrecht, The Netherlands

⁵INAF-Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone, Italy

⁶IRAP, Université de Toulouse, CNRS, CNES, 9 avenue du Colonel Roche, F-31028, Toulouse, France

⁷Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA

⁸Quasar Science Resources S.L for European Space Agency (ESA), ESAC, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain

⁹CRESST and X-ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

¹⁰Department of Physics and Astronomy, University of Southampton, Highfield, Southampton, SO17 1BJ, UK

¹¹Department of Physics and Astronomy, FI-20014, University of Turku, Finland

¹²Space Research Institute of the Russian Academy of Sciences, Profsoyuznaya Str. 84/32, Moscow 117997, Russia

¹³Space Sciences Laboratory, University of California, 7 Gauss Way, Berkeley, CA 94720-7450, USA

¹⁴Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

¹⁵Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield, AL10 9AB, UK

¹⁶Remeis-Observatory and Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Sternwartstr. 7, D-96049 Bamberg,
Germany

Received 2022 July 29; revised 2022 August 22; accepted 2022 August 26; published 2022 October 5

Abstract

M82 X-2 is the first pulsating ultraluminous X-ray source discovered. The luminosity of these extreme pulsars, if isotropic, implies an extreme mass transfer rate. An alternative is to assume a much lower mass transfer rate, but with an apparent luminosity boosted by geometrical beaming. Only an independent measurement of the mass transfer rate can help discriminate between these two scenarios. In this paper, we follow the orbit of the neutron star for 7 yr, measure the decay of the orbit ($P_{\text{orb}}/P_{\text{orb}} \approx -8 \cdot 10^{-6} \text{ yr}^{-1}$), and argue that this orbital decay is driven by extreme mass transfer of more than 150 times the mass transfer limit set by the Eddington luminosity. If this is true, the mass available to the accretor is more than enough to justify its luminosity, with no need for beaming. This also strongly favors models where the accretor is a highly magnetized neutron star.

Unified Astronomy Thesaurus concepts: Pulsars (1306); Ultraluminous X-ray sources (2164); Pulsar timing method (1305); Orbital evolution (1178)

1. Introduction

The luminosity of accreting sources is largely driven by the amount of matter that is transferred onto the accreting object, whether it be from a donor star for typical neutron stars and stellar-mass black holes, or an accretion disk for supermassive black holes at the centers of galaxies (Frank et al. 2002). There is a classical limit to the mass transfer, which corresponds to the mass-accretion rate that leads to a balance between the force of radiation pressure pushing outward and the gravitational force acting inward on an accreting object of mass *M*. For spherical hydrogen accretion, this corresponds to the Eddington luminosity:

$$L_{\text{Edd}} \approx 1.3 \cdot 10^{38} \frac{M}{M_{\odot}} \text{erg s}^{-1}. \quad (1)$$

Therefore, the extreme luminosity of ultraluminous X-ray sources (ULXs; Kaaret et al. 2017; Fabrika et al. 2021) led many to think that these sources were powered by intermediate-

mass black holes. Over the years, multiple pieces of evidence cast doubt on the applicability of this classical limit on ULXs (Poutanen et al. 2007; Gladstone et al. 2009; Bachetti et al. 2013). Eventually, the discovery of pulsating ultraluminous X-ray sources (PULXs; Bachetti et al. 2014, hereafter B14), accreting neutron stars radiating hundreds of times above their Eddington limits, demonstrated that super-Eddington accretion was a viable explanation for the majority of ULXs. It is still unclear how these pulsars (pulsating neutron stars) emit this extreme luminosity. Some argue that the isotropic luminosity is much lower, and the observed luminosity is boosted by geometrical beaming, driven by the collimation of a (less extreme) super-Eddington disk (King et al. 2017). This interpretation has found some support in global MHD simulations of accreting black holes and neutron stars, where mild-to-extreme geometrical beaming is observed (e.g., Jiang et al. 2014; Abarca et al. 2021). However, these simulations assume a low magnetic field of the neutron star ($\lesssim 10^{10}$ G), if any, and this collimation effect is likely to be lessened when the magnetic field of the pulsar is stronger. In fact, other models explain the luminosity with arguments centered on a high magnetic field of the pulsar ($>10^{13}$ G), like the reduction of the

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Thomson scattering cross section in high magnetic fields, either in their dipolar (Mushtukov et al. 2015, 2017) or their multipolar components (Brice et al. 2021). This reduction of the cross section allows to hit the local Eddington limit at much higher mass-accretion rates, increasing the maximum luminosity. It is also possible that the solution is a mixture of genuine super-Eddington accretion and a small amount of beaming (Israel et al. 2017).

A key difference between these models is the relation that they assume between the mass-accretion rate $\dot{m} = \dot{M}/\dot{M}_{\text{Edd}}$ and the luminosity, linear in the low-beaming scenario, almost quadratic ($L \propto (1 + \log \dot{m})\dot{m}^2$) in the other, due to the assumed quadratic dependence of beaming on the mass-accretion rate (King 2008). In other words, beaming models infer a much lower mass transfer rate between the donor star and the neutron star for a given luminosity.

An independent measurement of the mass transfer is key for disentangling these two scenarios. In principle, one way to measure this transfer of matter between two orbiting objects is through the observation of a decay of the orbital period (Tauris & van den Heuvel 2006).

One of the best systems where this can be tested is M82 X-2, the first PULX ever discovered. B14 and Bachetti et al. (2020, hereafter B20) measured the orbit of this PULX very precisely, determining an orbital period of 2.532948(4) days, a semimajor axis of 22.215(5) ls, and no detectable eccentricity ($<0.003$). What makes this system particularly interesting from the point of view of orbital decay measurements is that its revolution period is short enough, and the ephemeris known so precisely, that the epoch of passage through the ascending node can be constrained to $\sim 100$ s with a single, reasonably long X-ray observation.

In this paper, by tracking the ascending node passages over 8 yr of Nuclear Spectroscopic Telescope Array (NuSTAR) observations, we present the precise measurement of orbital decay in M82 X-2, leading to an estimate of mass transfer whose value agrees to within a factor 2 to the one inferred from the luminosity of the pulsar.

In Section 2 we describe the data reduction, and in Section 3 we detail the temporal analysis that led to the orbital decay measurement, while we devote the last sections to the interpretation of this orbital decay.

## 2. Data Reduction

### 2.1. NuSTAR

We downloaded all the NuSTAR data of M82 from the High Energy Astrophysics Science Archive Research Center (HEA-SARC). We ran *nupipeline* with standard options to produce cleaned event files. This tool produces different event files corresponding to different observing modes: *SCIENCE* (01), *OCCULTATION* (02), *SLEW* (03), *SAA* (04), *CALIBRATION* (05), and *SCIENCE_SC* (06). The modes usable for science are 01 and 06. Note, however, that only mode 01 data are recorded in normal instrumental conditions. Mode 06 data correspond to time intervals where only a subset of the camera head units (CHUs) are available, and the astrometry can be off by $1'-2'$ (see Walton et al. 2016 for an example of the astrometry issues in this observing mode). For mode 01 data, we used a region of 70" around the centroid of the X-ray source corresponding to the position of M82 X-1 and M82 X-2, which is spatially unresolved in NuSTAR. The centroid was

calculated independently for each observation and for each of the two focal plane modules, as a mismatch of $\sim 10''$ can be expected.

We processed mode 06 data with the *nusplitsc* tool, which separates events corresponding to different CHU combinations. For each of these event files, we adjusted the centroid of the source and repeated the selection done for mode 01 data. Finally, we merged the source-selected event lists from mode 01 and mode 06 data. In only a few cases, due to the source falling on a chip gap, we saw that the light curve showed visible “steps” between intervals corresponding to different CHU combinations. We verified that the addition or elimination of the problematic intervals did not alter significantly the power around the pulsation frequency $\sim 0.7$ Hz.

Finally, we ran *barycorr* to refer the photon arrival times to the solar system barycenter. We selected the ICRS reference frame, the DE421 JPL ephemeris, and the position of M82 X-2 determined by Chandra. For all observations, we used the latest clock-correction file available, which provides an absolute time precision of $\sim 60$ µs.

## 3. Timing Analysis

Due to the very low pulsed fraction in the XMM-Newton band, demonstrated in Appendix F, we only used NuSTAR data for timing analysis.

Initially, we largely followed the search strategies used in B20, running $Z_I^2$ searches (Buccheri et al. 1983), also known as the Rayleigh test, on the event arrival times corrected for orbital motion, varying the ascending node passage epoch $T_{\text{asc}}$ on a fine grid between $-P_{\text{orb}}/5$ and $P_{\text{orb}}/5$. This time, the search allowed a range of spin derivatives for each trial ascending node passage value. The spin parameters vary so rapidly that they are only loosely constrained when observations are just $\sim 1$ week apart. There is no way to reliably phase-connect separate observations. Therefore, even observations done $\sim 2$ weeks apart were analyzed singularly. For the search in the *f-f* plane, we used the “quasi-fast folding algorithm” (B20), which calculates the $Z_I^2$ on prebinned profiles (Bachetti et al. 2021), using at least 16 bins for the folded profiles. Moreover, we ran the search both in the full-energy band and between 8 and 30 keV. This allowed four detections in the new observations. We also reran a pulsation search using all available observations, and we found highly significant ($\ge 5\sigma$) pulsations in two archival data sets, corresponding to ObsIDs 30101045002 and 80002092002. In both observations, pulsations are more strongly detected in the 8–30 keV energy band. Moreover, surprisingly, we find that during 30101045002 the pulsar was instantaneously spinning *down*. This is the first time that M82 X-2 is found to be spinning down while accreting and pulsating, and provides clear evidence that a significant part of the torque from the disk is happening outside the corotation radius (see Appendix C for more details). This is probably why B20 only obtained marginal evidence for pulsations in this ObsID. This new detection is important because pulsations are detected over a $\sim 4$ day interval, which is long enough to provide an excellent constraint on $T_{\text{asc}}$. For each detection, we then ran the search again around the best solution, oversampling by a large factor to find the best estimate of the mean of each parameter.

We proceeded to create local timing solutions for each observation, leaving the orbital parameters from B20 unchanged with the exception of the ascending node, and
<table><caption>**Table 1**
Spin and Orbital Parameters from the Multidimensional Timing Procedure in Section 3</caption><thead><tr><th>Obs. ID</th><th>Epoch (MJD)</th><th>Energy (keV)</th><th>T<sub>asc</sub> (MJD)</th><th>f<sub>spin</sub> (Hz)</th><th>f&#x0307;<sub>spin</sub> (10<sup>-12</sup> Hz s<sup>-1</sup>)</th><th>f&#x0307;&#x0307;<sub>spin</sub> (10<sup>-16</sup> Hz s<sup>-2</sup>)</th><th>&Delta;T<sub>asc</sub> (s)</th></tr></thead><tbody><tr><td>80002092002*</td><td>56681.24441</td><td>8–30</td><td>56682.073(5)</td><td>0.728509(4)</td><td>-90(70)</td><td></td><td>600(500)</td></tr><tr><td>80002092004</td><td>56683.81009</td><td>8–30</td><td>56684.6004(28)</td><td>0.7285316(16)</td><td>50(35)</td><td></td><td>110(240)</td></tr><tr><td>80002092006</td><td>56688.80899</td><td>8–30</td><td>56689.6656(5)</td><td>0.72854791(22)</td><td>15.3(12)</td><td>1.20(27)</td><td>50(40)</td></tr><tr><td>80002092007</td><td>56694.12259</td><td>3–30</td><td>56694.73184(19)</td><td>0.72856174(6)</td><td>34.6(14)</td><td></td><td>86(17)</td></tr><tr><td>80002092007</td><td>56697.38070</td><td>3–30</td><td>56697.2648(4)</td><td>0.72857943(8)</td><td>92(6)</td><td></td><td>90(40)</td></tr><tr><td>80002092008</td><td>56700.75316</td><td>3–30</td><td>56699.798(5)</td><td>0.728609(5)</td><td>90(60)</td><td></td><td>100(400)</td></tr><tr><td>80002092009</td><td>56700.75316</td><td>3–30</td><td>56699.7978(4)</td><td>0.72860925(13)</td><td>106(5)</td><td></td><td>90(32)</td></tr><tr><td>80002092011</td><td>56720.87754</td><td>8–30</td><td>56720.0612(12)</td><td>0.7287596(6)</td><td>113(13)</td><td></td><td>80(100)</td></tr><tr><td>30101045002*</td><td>57495.31178</td><td>8–30</td><td>57495.1440(7)</td><td>0.72519103(20)</td><td><b>-65(5)</b></td><td></td><td>140(60)</td></tr><tr><td>90201037002</td><td>57641.99852</td><td>8–30</td><td>57642.049(8)</td><td>0.7239040(12)</td><td>-320(230)</td><td></td><td>-400(700)</td></tr><tr><td>30502021002*</td><td>58919.09530</td><td>3–30</td><td>58918.6261(12)</td><td>0.7219294(7)</td><td>36(15)</td><td></td><td>-2860(100)</td></tr><tr><td>30602027002*</td><td>59312.65089</td><td>3–30</td><td>59313.750(6)</td><td>0.7222096(21)</td><td>120(150)</td><td></td><td>-4300(500)</td></tr><tr><td>30602027004*</td><td>59326.05680</td><td>8–30</td><td>59326.409(4)</td><td>0.7222978(25)</td><td>50(70)</td><td></td><td>-4700(400)</td></tr><tr><td>30702012002*</td><td>59505.27806</td><td>3–30</td><td>59506.2428(11)</td><td>0.72086594(33)</td><td><b>-45(14)</b></td><td></td><td>-5210(90)</td></tr></tbody></table>

Note. Starred ObsIDs are those corresponding to new detections from this paper. We also highlight in bold observations with significant (>3σ) evidence of spin down. The energy range is the one where the pulsations are detected with the highest significance. Data from ObsID 80002092006 start after the glitch reported by B20 at MJD 56685.7. ObsID 80002092007 has a sudden change of frequency, probably another glitch, around MJD 56696, therefore we split the observation in two parts around that epoch. ΔT<sub>asc</sub> was calculated with respect to the orbital ephemeris from B20.

setting the current spin frequency and frequency derivative. These local solutions differed only for the parameters F0 (spin frequency in Hertz), F1 (spin derivative in Hertz per second), and TASC (epoch of passage through the ascending node in Modified Julian Date, MJD). Then, we used the method by Pletsch & Clark (2015) to make a Bayesian fit of the local timing solution, using a sinusoidal pulse template normalized to the same pulsed fraction of the pulsar in each given observation. Working in phase space instead of frequency, this method is far more sensitive to small changes of parameters, and yields very precise estimates on them. The exact parameters we fitted were the difference from F0 in units of 10<sup>-x</sup> s, depending on the observation length, the difference from F1 in units of 10<sup>-y</sup> Hz s<sup>-1</sup> (where X, Y were chosen as values close to the order of magnitude of the known errors on the parameters), and the difference from TASC in seconds. This was done to avoid incurring any numerical errors due to the small steps involved in some parameters. We set flat priors for all parameters: 0.5 < f < 1 Hz, |f| < 10<sup>-7</sup>, and ΔT<sub>asc</sub> < P<sub>orb</sub>. We first used the scipy.optimze.minimize function to minimize the negative log-likelihood and determine an approximate starting solution. Then, values around this solution were used to initialize a Markov Chain Monte Carlo (MCMC) sampler as implemented in the emcee (Foreman-Mackey et al. 2013) library. Since the analysis took a significant time (up to 2 s per iteration in the larger data sets), we followed the instructions in the emcee documentation to interrupt the sampling once the iterations had reached 200 times the “autocorrelation time” τ, more than the recommended 50 for additional robustness. τ itself was calculated every 100 steps of the chain. The number of steps used in the chains varied between 3000 and 100,000 depending mostly on the length of the observation and the number of photons, with longer observations requiring fewer steps (because of the reduced correlations between parameters). We used the 3–30 keV or 8–30 keV energy range depending on which range yielded the highest power in the Rayleigh search. This allowed us to estimate the posterior distribution on the parameters and their uncertainties. The posterior distributions

are generally well behaved, with reasonably (sometimes slightly skewed due to the correlation between the parameters) bell-shaped distributions. We determined 1σ error bars on the parameters by looking at the 16% and 84% percentiles. The results are summarized in Table 1, and the detection in ObsID 30101045002 is shown in Figure 1 as an example.

### 3.1. Orbital Decay

To describe a change of the orbital period over time, it is customary to measure the time that a star reaches a particular phase of the orbit, and compare it with the expected time given the previous orbital solution. For circular orbits with no eclipses, it is common to use one of the two intersections between the orbit and the plane perpendicular to the line of sight passing through the center of mass of the binary system. These points of the orbit are called nodes; the ascending node is the node that the pulsar crosses when moving away from the observer. The expected time of passage at the ascending node after *n* orbits, T<sub>asc,n</sub> (for simplicity, we drop the *asc* when *n* or other indices are present), in the presence of an orbital derivative, can be expressed as (see Kelley et al. 1980; Falanga et al. 2015)

$$T_n = T_0 + n P_{\text{orb}} + \frac{1}{2} n^2 P_{\text{orb}} \dot{P}_{\text{orb}} + \frac{1}{6} n^3 P_{\text{orb}}^2 \ddot{P}_{\text{orb}} + \dots \quad (2)$$

By using a previously determined ephemeris as a baseline, we can measure the delay of the measured T<sub>asc</sub> from the expected one. When plotting this delay, offsets indicate a shift in T<sub>asc</sub>, linear trends an uncertainty of P<sub>orb</sub>, and parabolic trends an orbital period derivative:

$$\delta T_n(t) = \delta T_{\text{asc}} + \frac{t - T_{\text{asc}}}{P_{\text{orb}}} \delta P_{\text{orb}} + \frac{1}{2} \frac{\dot{P}_{\text{orb}}}{P_{\text{orb}}} (t - T_{\text{asc}})^2, \quad (3)$$

where we substituted $n = (t - T_{\text{asc}})/P_{\text{orb}}$.

Using the new T<sub>asc</sub> values, we infer the orbital decay of M82 X-2 using a Bayesian model.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Example detection of pulsations, and orbital/spin parameter refinement, from ObsID 30101045002. The color map shows the Rayleigh search in a three-parameter grid, using spin frequency, spin first derivative, and the drift of the periastron passage $T_{asc}$. The corner plot on the right, instead, shows the refinement of the results of the Rayleigh search, with the addition of pulse phase, fitted using the Bayesian method by Pletsch & Clark (2015). For example, if $dF0-5$ is -0.03 and the initial F0 was 0.725, this means that the best-fit frequency is $0.725 - 0.03 \cdot 10^{-5}$ Hz.</figcaption></figure>

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

The following equation serves as the orbital evolution model:

$$ \Delta T_{asc} = a + \frac{b}{P_{orb}}(t - T_{asc}) + 0.5 \cdot 10^{-6} c \frac{86400}{365.25}(t - T_{asc})^2, \quad (4) $$

where $t$, $T_{asc}$, and $P_{orb}$ are expressed in days, $\Delta T_{asc}$ is the delay of $T_{asc}$ in seconds, $a$ is a correction to $T_{asc}$ in seconds, $b$ is a correction to $P_{orb}$ in seconds, and $c$ is the new value of $\dot{P}_{orb}/P_{orb}$ in units of $10^{-6}$yr$^{-1}$. The baseline solution from B20 was $T_{asc,B20} = \text{MJD} 56682.0661$, $P_{orb,B20} = 2.532948$ days, and $\dot{P}_{orb} = 0$ s s$^{-1}$. Priors for $a$, $b$, and $c$ were uniform between $\pm 10^6$; in checks, we found that the width of the prior has no significant effect on our posterior inference.

We first performed a maximum a posteriori fit with a standard Gaussian likelihood, allowing for asymmetric error bars. The solution served as an initialization of a MCMC sampler using `emcee`, as before. Using 32 walkers, we ran the chains for 20,000 steps. We calculated the autocorrelation “time,” which was at most 46 steps. We thinned the chain by a factor 23 (half the autocorrelation length) and discarded 920 steps (20 times the autocorrelation length) as burn-in. The resulting marginal posterior probability distributions are plotted using the `corner` library (Foreman-Mackey 2016) in Figure 2.

We find posterior means and credible intervals of $a = 72$ (13), $b = 2.18(26)$, and $c = -8.20(34)$. Using these values, we corrected the orbital parameters as $T_{asc} = T_{0,B20} + a$ sec, $P_{orb} = P_{orb,B20} + b$ sec, and $\dot{P}_{orb}/P_{orb} = c \cdot 10^{-6}$yr$^{-1}$ to obtain the values in Table 2.

Finally, we fixed the orbital parameters and we reran a final accelerated search for pulsations in all ObsIDs using the Rayleigh test, yielding the results in Table 1.

# 4. Discussion

Over the years, many models have been proposed to describe the interaction between the plasma in the disk and the magnetic field lines of an accreting pulsar (Ghosh & Lamb 1978; Wang 1996; Chashkina et al. 2017). Despite large differences in the treatment of the details of this interaction, these models make estimates for the magnetic field within ~one order of magnitude when the inner radius and the mass-accretion rate are fixed (Xu & Li 2017; Erkut et al. 2020; Chen et al. 2021), if one can assume spin equilibrium: a regime where the outward pressure from the rotating magnetic field balances almost exactly the ram pressure from the infalling matter.

Until now, different groups have used the observed luminosity as a proxy for the mass-accretion rate, and this produced very different estimates depending on the assumption of the beaming fraction. In addition, different works used different assumptions on the position of the inner radius, with the high-magnetic-field models assuming spin equilibrium (Dall’Osso et al. 2015; Eksi et al. 2015; Tsygankov et al. 2016) and the beaming models being incompatible with it (King et al. 2017).

In this work, we produce robust evidence in favor of spin equilibrium, and we measure an orbital decay that might provide an independent estimate of mass transfer, as we are going to discuss below.

## 4.1. Spin Equilibrium

Thanks to the new detections listed above, we found that for at least part of the time between 2016 and 2020 the pulsar continued to spin down (slow down its rotation), as reported by B20, because the frequency (~0.721 Hz) observed in 2020 was lower than observed in 2016 (~0.723 Hz). However, since then, the neutron star appears to be alternating phases of spin up and spin down around ~0.721 Hz. In at least one observation in 2016 and probably in another in 2021, the
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Left: orbital decay in M82 X-2, measured through the delay of the $T_{asc}$ parameter (time of passage through the ascending node) from Equation (3). Semitransparent orange lines indicate possible ~1000 quadratic solutions, coming from MCMC sampling (described in Section 3.1). Grey points are lower-quality measurements (uncertainty in $T_{asc}$ larger than 200 s due to high orbital-spin parameter correlation). Right: corner plot of the posterior distribution of the parameters of the orbital decay with Equation (4), sampled with MCMC. Vertical dashed lines show the 16%, 50%, and 84% percentiles.</figcaption></figure>

<table><caption>**Table 2**
Updated Orbital Parameters for M82 X-2, as Determined in This Work</caption><thead><tr><td>Parameter</td><td>Unit</td><td>Value (Uncert)</td></tr></thead><tbody><tr><td>P<sub>orb</sub></td><td>d</td><td>2.5329733(30)</td></tr><tr><td>P&#x0302;<sub>orb</sub></td><td>s s<sup>-1</sup></td><td>-5.69(24) &middot; 10<sup>-8</sup></td></tr><tr><td>P&#x0302;<sub>orb</sub>/P<sub>orb</sub></td><td>yr<sup>-1</sup></td><td>-8.20(34) &middot; 10<sup>-6</sup></td></tr><tr><td>a sin i</td><td>light-seconds</td><td>22.218(5)</td></tr><tr><td>T<sub>asc</sub></td><td>MJD</td><td>56682.06694(15)</td></tr><tr><td>e</td><td></td><td>&lt;0.0015 (3σ u.l.)</td></tr></tbody></table>

pulsar was spinning down while accreting (see Table 1). In summary, the spin evolution of M82 X-2 strongly points to a situation of spin equilibrium. In this condition, spin up and spin down can be produced with small changes of accretion rate (D’Angelo & Spruit 2012), and it is possible to confidently estimate the magnetic field of the neutron star by equating the analytical formulas for the inner radius $R_{in}$ to the corotation radius $R_{co}$, at which the angular velocity of the matter in the disk equals the one of the star (see Appendix C).

Being close to spin equilibrium also implies that a relatively small drop of mass-accretion rate could trigger the so-called “propeller” regime (Illarionov & Sunyaev 1975), where the rotating, highly magnetized pulsar is able to swipe away the infalling matter. During the transition to this regime, it is possible to still have accretion, albeit discontinuous (Romanova et al. 2004), up to a point where accretion is stopped altogether, leaving only the disk and the ejected matter as sources of X-ray flux. Based on a possible bimodal distribution of the fluxes of M82 X-2, Tsygankov et al. (2016) claimed that the observed low states in M82 X-2 were evidence of this transition. It is not clear, at the moment, if this is compatible with the observed ~60 day periodicity of the low states (Brightman et al. 2019), which would imply a *periodic* decrease of mass transfer, difficult to reconcile with the very low eccentricity of the system.

## 4.2. Is it Mass Transfer?

The observed orbital decay is compatible with the mass transfer from a more massive donor star to a neutron star (Tauris & van den Heuvel 2006; see Appendix B). Assuming a pulsar mass $M_p = 1.4 M_\odot$ and a donor mass $M_d = 8 M_\odot$ (which corresponds to the mean of the probability distributions of masses; see Appendix D), it is straightforward to estimate the mass transfer rate from the observed orbital decay, assuming conservative mass transfer, as $\dot{M}_d \approx -4.7 \cdot 10^{-6} M_\odot \text{ yr}^{-1}$. This corresponds to ~200 times the Eddington limit, assuming an Eddington mass accretion rate corresponding to $\dot{M}_{Edd} \approx 1.5 \cdot 10^{18} \text{ g s}^{-1} \approx 2.4 \cdot 10^{-8} M_\odot \text{ yr}^{-1}$. This is the mass that the donor transfers into the Roche lobe of the neutron star. This mass exchange exceeds both that inferred from the apparent bolometric luminosity of the source (which is at most ~100 times the Eddington limit; see B14), or the one inferred from beaming scenarios (36 times Eddington; King et al. 2017); see Figure 3 for details. It is possible that part of this matter leaves the system through fast winds launched from the super-Eddington disk (Pinto et al. 2016; Kosec et al. 2018). This mass loss happens from the vicinity of the accretor, and its specific angular momentum is such that the effect on the orbit is similar to the conservative case (see Appendix B). On the other hand, it decreases the amount of matter that accretes onto the neutron star, which is a viable explanation for the slightly lower luminosity observed. Isotropic mass loss from the donor, instead, as one would expect from stellar winds, would have the opposite effect, expanding the orbit. Additional mass loss from the outer disk in the form of slow winds (e.g., Middleton et al. 2022) would represent an intermediate case, carrying away specific angular momentum somewhere between that at the position of the neutron star and that at the first Lagrangian point, L1. Therefore, the estimate above is a *lower limit* to the mass transfer rate. Another possibility, involving mass loss from the second Lagrangian point, L2, forming a circumbinary disk, is discussed below.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Mass transfer rate toward the accretor (through L1), vs. donor mass, given the measured orbital period decay. We estimate the mass transfer rate for different mechanisms of mass loss and different accretor masses. Solid lines indicate conservative mass transfer, dashed lines indicate 90% isotropic mass loss from the donor, and dotted lines indicate 100% mass loss from the proximity of the accretor. The latter is only a small correction to the conservative case, while mass loss from the donor implies a much larger mass transfer rate in order to produce the observed orbital decay. Vertical bands on the left show the limit donor masses given the absence of eclipses, for different accretor masses (same color coding).</figcaption></figure>

If the observed orbital decay really is due to mass transfer, we can fix these two important variables, leading to an interpretation of M82 X-2 being a highly magnetized neutron star ($B > 10^{13}$ G) with any of these models (see also Appendix E), and exclude beaming as a primary amplifier of ULX emission (see also Vasilopoulos et al. 2021 for further evidence in this sense).

## 4.3. Alternative Models

### 4.3.1. Synchronization and Circularization

The value of orbital derivative found in this work is larger, albeit only by a few times, than those observed in much less luminous X-ray pulsars such as SMC X-1, LMC X-4, and Cen X-3 (-3.5, -1.0, and -1.8 $\cdot 10^{-6}$ yr$^{-1}$, respectively; see Falanga et al. 2015). The slight mismatch that we find in M82 X-2 between the mass transfer inferred from the orbit and that inferred from the luminosity becomes a factor ~20 for SMC X-1. This has led Falanga and many other authors (e.g., Levine et al. 1993, 2000; Chernov 2020) to disregard mass transfer as the primary engine for orbital decay in these systems. At this point, we cannot exclude that whatever process believed to be in place in those systems (like viscous processes producing the circularization of an elliptical orbit or the synchronization of the star's rotation with the orbit; e.g., Falanga et al. 2015, and references therein; Chernov 2020) is at work in M82 X-2. However, we do stress that these high-mass X-ray binaries (HMXBs) are likely accreting in very different regimes, possibly from focused winds and not from a Roche-lobe overflow, as would instead be expected from a super-Eddington source.

### 4.3.2. Circumbinary Disk

The observed orbital decay is, in principle, compatible with an equatorial circumbinary disk launched by L2 (Tauris & van den Heuvel 2006; see also Appendix B). However, this only

happens in situations where the donor star inflates well beyond its Roche lobe, or via a fast and stable wind, and its onset quickly leads to unstable orbital decay and common envelope (Misra et al. 2020). Lu et al. (2022) studied in detail the conditions for this phenomenon, finding that it should be important above $\dot{m} \sim 10^{-4} M_{\odot} \text{yr}^{-1}$. SS433, a possible ULX analog in our Galaxy (Fabrika et al. 2006; Middleton et al. 2021), might be undergoing such a process. However, in that case, the accretor is believed to be a stellar-mass black hole (Blundell et al. 2008) and the mass ratio is ~1, and this process leads to the expansion of the orbit (Cherepashchuk et al. 2021), stabilizing the mass transfer.

## 5. Conclusion

The detection of orbital decay in M82 X-2 is a key milestone to understand the evolution of this system and, possibly, of all low-orbital-period PULXs like NGC 5907 ULX1 (Israel et al. 2017) and M51 ULX-7 (Rodríguez Castillo et al. 2020).

We argue that the decay is driven by mass transfer: the implied mass transfer is only a factor ~2 above the one inferred from the luminosity, and this can easily be explained by a slightly lower efficiency or a massive outflow such as those observed in other ULXs (but undetectable in M82 X-2 due to source confusion in the M82 field). Currently, we cannot exclude that phenomena such as the synchronization of the donor rotation with the orbit and/or the circularization of the orbit are contributing to the observed decay, which is only a few times higher than that of the eclipsing HMXBs from Falanga et al. (2015). Note, however, that this source has a quite tight upper limit on the eccentricity (see Appendix A) and the accretion is very likely through Roche-lobe overflow, at a much higher rate than the sample from Falanga et al. (2015), making the timescale for synchronization faster. If these phenomena are in place, it is likely that we are witnessing a very-short-lived phase of the evolution of this binary system.

Regardless of the exact driver of the observed orbital decay, our measurement informs the theoretical study of ULX progenitors. At the moment, the evolutionary scenarios able to produce a ULX seem often to lead to common envelope, with a relatively short phase of extreme mass transfer, in particular for donor masses in the lower mass of the allowed ranges for M82 X-2 (Tauris & van den Heuvel 2006; Misra et al. 2020). Stabilizing mechanisms, such as mass loss from the donor, are often invoked to increase the lifespan of ULXs, provided that the envelope is radiative.

We encourage further theoretical studies on the evolution of binary systems, to understand the conditions in which an orbital decay such as the one we observe can be produced. Future missions with instruments at higher throughput, like Athena (Barcons et al. 2017), will help detect pulsations and perform similar studies in many more ULXs. For M82 X-2 and extragalactic pulsars in general, which have hard pulsations and are often found in crowded fields, hard imagers with high angular resolution and good timing capabilities, like the proposed NASA probe High-Energy X-ray Probe (Madsen et al. 2019), would be excellent. Timing-devoted missions with large collecting area, such as the Chinese-Italian enhanced X-Ray Timing and Polarimetry mission (Zhang et al. 2016) or the proposed NASA probe Spectroscopic Time-Resolving Observatory for Broadband Energy X-rays (Ray et al. 2018), will allow sensitive searches for pulsations and timing studies in ULXs, provided that they are sufficiently isolated.
The authors wish to thank Victoria Grinberg, Włodek Kluźniak, and Alessandro Ridolfi for useful discussions, and the staff at the NuSTAR Science Operations Center at Caltech for the help in scheduling the observations and the frequent clock-correction file updates, which allowed a prompt analysis of the data. We would also wish to thank the anonymous referee, and the three referees of a previous submission, who provided very insightful feedback that led to a substantial improvement of the quality of the analysis. M.B. was funded in part by PRIN TEC INAF 2019 “SpecTemPolar!—Timing analysis in the era of high-throughput photon detectors”. M.H. is supported by an ESO fellowship. G.L.I. and M.B. acknowledge funding from the Italian MIUR PRIN grant No. 2017LJ39LM. A.D.J. was funded in part by the Chandra grant No. 803-0000-716015-404H00-6100-2723-4210-40716015HH 83121. J.P. was supported by the grant No. 14.W03.31.0021 of the Ministry of Science and Higher Education of the Russian Federation and the Academy of Finland grant No. 333112. D.J. W. acknowledges support from STFC in the form of an Ernest Rutherford Fellowship. H.P.E. acknowledges support under NASA Contract No. NNG08FD60C.

All the analysis of this paper was done using open-source software, Astropy, Stingray, HENDRICS, PINT, emcee, corner, and scinum, and can easily be verified using the solutions in Tables 1 and 2. The implementation of the Pletsch & Clark (2015) method can be found in the github repository https://github.com/matteobachetti/ell1fit. Figures were produced using the Matplotlib library and the Veusz software. The data used for this work come from the NuSTAR and XMM-Newton missions and are usually held private for one year, and made public on the High Energy Astrophysics Science Archive (HEASARC) and the XMM-Newton Science Archive (XSA) afterwards. NuSTAR is a Small Explorer mission led by Caltech and managed by the JPL for NASA's Science Mission Directorate in Washington. NuSTAR was developed in partnership with the Danish Technical University and the Italian Space Agency (ASI). The spacecraft was built by Orbital Sciences Corp., Dulles, Virginia. XMM-Newton is an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.

*Facilities:* NuSTAR, XMM.

*Software:* astropy (Astropy Collaboration et al. 2018), Stingray (Huppenkothen et al. 2016, 2019), HENDRICS (Bachetti 2018), PINT (Luo et al. 2021), scipy (Virtanen et al. 2020), numpy (Harris et al. 2020), numba (Lam et al. 2015), Veusz,<fnref n="17" /> Matplotlib (Hunter 2007), corner (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013).

## Appendix A
Eccentricity and Semimajor Axis

Thanks to the additional counts coming from the data reduction described in Section 2, we could obtain a more stringent upper limit on the eccentricity and verify the past estimates of the semimajor axis.

We created a piecewise spindown solution for PINT (using the PiecewiseSpindown model) using all the best estimates of the frequency and the frequency derivative listed in Table 2, which served as a baseline for the subsequent calculations. We used HENphaseogram (Bachetti 2018) to obtain times of arrival (TOAs) in 120 high signal-to-noise time intervals between MJD 56685.7 and 56722. Then, we used

pintk to look for features in the timing residuals reminiscent of an eccentricity. The Roemer delay gives the delay of the signal from the pulsar during its binary motion. In the limit of small eccentricity, this delay can be expressed as

$$ \Delta_R = X \left[ \sin \Phi + \left( \frac{\kappa}{2} \sin 2\Phi + \frac{\eta}{2} \cos 2\Phi \right) \right], \quad (A1) $$

where $X = a \sin i/c$, $\Phi = 2\pi/P_{\text{orb}}(T - T_{\text{asc}})$, $\omega$ is the angle of periastron, $\kappa = e \sin \omega$, and $\eta = e \cos \omega$. Hence, eccentricity should produce sinusoidal residuals with $P = P_{\text{orb}}/2$ and amplitude $e a \sin i/c$ in the TOAs corrected with a circular orbit. These features are not correlated with any other orbital parameter of interest (which produce features at the orbital period), and can be investigated independently. The Tempo2/PINT timing model ELL1 (Lange et al. 2001) implements this correction. Using PINT, we fit the best-fit residual from the best circular model with ELL1 and found no significant features reminiscent of an eccentricity. The new $3\sigma$ upper limit on eccentricity, using the 0.0005 error bars from this fit, is around 0.0015, half the value quoted by B14.

Using (A1), it is also possible to compare the effect of an error on $a \sin i/c$ with that on $T_{\text{asc}}$. Neglecting the eccentricity, we get that a given error on $\Delta_R$ can be written as

$$ \delta\Delta_R \approx \delta X \sin \Phi + \delta T_{\text{asc}} \frac{2\pi X}{P_{\text{orb}}} \cos \Phi. \quad (A2) $$

For M82 X-2, $a \sin i/c = 22.218(5)$. Therefore, we can neglect the error on $X$ whenever the error on $T_{\text{asc}}$ satisfies

$$ \delta T_{\text{asc}} \gg \frac{\delta X}{X} \frac{P_{\text{orb}}}{2\pi} \approx 8 \text{ s.} \quad (A3) $$

This is always true in this paper (see Table 1). Later observations are not able to constrain both $T_{\text{asc}}$ and $X$, and thawing $X$ in the fit artificially increases the error bars without leading to a more precise estimate: for short observations, it correlates with $T_{\text{asc}}$ and $\dot{v}$ and the fit yields unreasonable values both for $X$ and the other parameters.

## Appendix B
Mass Transfer

By differentiating the formula for the orbital angular momentum and Kepler's third law, it can be shown how the orbital separation and the orbital period change as a response to mass transfer or angular momentum changes (e.g., Tauris & van den Heuvel 2006):

$$ \frac{2}{3} \frac{\dot{P}_{\text{orb}}}{P_{\text{orb}}} = \frac{\dot{a}}{a} = 2 \frac{\dot{J}}{J} - 2 \frac{\dot{M}_d}{M_d} - 2 \frac{\dot{M}_p}{M_p} + \frac{\dot{M}_d + \dot{M}_p}{M_d + M_p} - 2 \dot{e} e, \quad (B1) $$

where $J$ is the total angular momentum of the system, $M_p$ and $M_d$ are the masses of the pulsar and the donor, $a$ is the orbital separation, $e$ is the eccentricity, and dots denote time derivatives. $\dot{M}_d$ is negative and is $\dot{M}_p$ positive, because the pulsar is accreting from the donor. The eccentricity of M82 X-2 is consistent with 0 (see Appendix A), as expected from a Roche-lobe-overflowing system, so it is likely that the last term in the equation can be neglected. But admitting that an undetected tiny eccentricity exists in the system, in order to have a negative orbital derivative there should be a positive

<footnote n="17">¹⁷ https://veusz.github.io/</footnote>
change, or an increase, of eccentricity, which is implausible given that these systems tend to circularize over time.

A number of phenomena causing changes in orbital angular momentum are discussed in the literature, such as gravitational wave (GW) emission (important in very compact systems such as some binary neutron stars), spin-orbit coupling (when the Roche-filling star's rotation is not synchronized with the orbit), magnetic braking (studied in low-mass X-ray binaries), and mass loss when the ejected mass has specific angular momentum. Given the large donor mass and orbital distance, we do not expect GW emission or magnetic braking to be significant. Moreover, even though they disagree on the exact mass transfer rate, different authors agree that the system is undergoing a strong mass transfer (Mushtukov et al. 2017; King & Lasota 2020). Such a mass transfer rate is difficult to reconcile with mechanisms other than Roche-lobe overflow (such as wind accretion or even wind Roche-lobe overflow; Mellah et al. 2019), and the synchronization timescales are so small that we can also neglect spin-orbit coupling (Tauris & Savonije 2001; Stoyanov & Zamanov 2009). This leaves us with mass transfer and or mass loss from a circumbinary disk (see below) as the only likely sources of angular momentum drain.

Conservative mass transfer has no angular momentum or mass losses from the system (i.e., $\dot{M}_p = -\dot{M}_d$ and $\dot{J} = 0$). In this case, Equation (B1) reduces to

$$ \frac{\dot{P}_{\text{orb}}}{P_{\text{orb}}} = 3 \frac{\dot{M}_d}{M_d} \left( \frac{M_d}{M_p} - 1 \right). \qquad (B2) $$

It is clear that, for $M_d/M_p > 1$, the system responds to a mass transfer from the donor ($\dot{M}_d < 0$) by decreasing the orbital period, as observed.

The nonconservative mass transfer case (when mass is lost from the system in any form) implies a change of the total angular momentum and can be studied by dividing the angular momentum term into different terms. Following the approach by van den Heuvel (1994), Soberman et al. (1997), and Tauris & van den Heuvel (2006),

$$ \frac{\dot{J}}{J} = \frac{\alpha + \beta r^2 + \delta \gamma (1 + r)^2 \dot{M}_d}{1 + r} \frac{\dot{M}_d}{M_d} \qquad (B3) $$

and

$$ \dot{M}_p = -(1 - \alpha - \beta - \delta) \dot{M}_d, \qquad (B4) $$

where $r = M_d/M_p$, $\alpha$ indicates the fraction of matter lost directly from the donor, $\beta$ the fraction lost from fast winds close to the accretor, and $\delta$ the fraction lost in a circumbinary disk of radius $a_r = \gamma^2 a$.

It is interesting to show where the three angular momentum losses lead when they dominate the orbital evolution, by developing Equation (B1) with Equations (B3) and (B4).

For the loss from the donor ($\alpha = 1$):

$$ \frac{\dot{P}_{\text{orb}}}{P_{\text{orb}}} = \frac{3}{2} \frac{\dot{M}_d}{M_d} \left( \frac{-r}{1+r} \right). \qquad (B5) $$

Therefore, an isotropic mass loss from the donor leads to an expansion of the orbit.

For the loss from the accretor ($\beta = 1$):

$$ \frac{\dot{P}_{\text{orb}}}{P_{\text{orb}}} = \frac{3}{2} \frac{\dot{M}_d}{M_d} \left( \frac{2r^2 - r - 2}{1 + r} \right), \qquad (B6) $$

implying that isotropic mass loss from the accretor (e.g., with disk winds) still leads to a contraction of the orbit. This is what is believed to happen at extreme mass transfer rates, where we expect strong radiation-driven winds to be launched inside the spherization radius (Shakura & Sunyaev 1973). In the limit $r \gg 1$, this is equivalent to the conservative case.

Finally, for the circumbinary disk ($\delta = 1$), we have

$$ \frac{\dot{P}_{\text{orb}}}{P_{\text{orb}}} = \frac{3}{2} \frac{\dot{M}_d}{M_d} \left( \frac{2\gamma(1+r)^2 - 2 - r}{1+r} \right), \qquad (B7) $$

which, for $\gamma \gg 1$ (disk radius larger than orbital separation) and $r > 1$, also produces a contraction of the orbit.

To summarize, the orbital decay we observe is compatible with the effect of mass transfer between a more massive donor and a neutron star (with or without mass loss from the accretor), or with angular momentum loss through an equatorial circumbinary disk, possibly launched by L2 (Tauris & van den Heuvel 2006). Due to the observation that matter is indeed accreting onto the neutron star, and that we observe many ULXs in nearby galaxies which suggests that this accretion regime is not too short lived, our analysis favors conservative (or mildly nonconservative) mass transfer from an intermediate-/high-mass star, with no high-angular-momentum mass-loss mechanisms. Again, we stress that fast winds from the region around the compact object do not change the results considerably. Moreover, as we show in Appendix E, the spherization radius is likely in the proximity of the magnetospheric radius, changing these estimates by a relatively small amount.

## Appendix C: Important Radii

Around an accreting neutron star, we can define two important radii (see Frank et al. 2002 for a comprehensive treatment): the first, the corotation radius $R_{\text{co}}$, is the radius at which the Keplerian angular velocity in the disk equals the angular velocity of the neutron star:

$$ R_{\text{co}} = \left( \frac{G M_p \dot{M}_p^2}{4 \pi^2} \right)^{\frac{1}{3}}, \qquad (C1) $$

where $p_{\text{spin}}$ is the rotation period of the neutron star, $M$ its mass, and $G$ the universal gravitational constant.

The second is called the magnetospheric radius, or inner radius, or truncation radius. Within this radius, the accretion disk gets disrupted, and matter gets captured by the magnetic field lines and conveyed to the magnetic poles of the neutron star:

$$ R_{\text{in}} = \xi \left( \frac{\mu^4}{G M \dot{M}^2} \right)^{\frac{1}{7}}, \qquad (C2) $$

where $\mu$ is the magnetic dipole moment, $\dot{M}$ the mass-accretion rate, and $\xi \sim 0.5$ encodes a number of effects like the accretion geometry (e.g., disk versus isotropic accretion) and the details

<footnote n="18">¹⁸ Note that in other papers (e.g., Joss & Rappaport 1984), $\alpha$ indicates the specific angular momentum. This can create confusion when comparing the different approaches.</footnote>
of the interaction between the plasma and the different components of the magnetic field.

According to accretion theory, the relative position of $R_{co}$ and $R_{in}$ is what determines whether a neutron star will spin up (accelerate its rotation) during accretion or spin down (slow down its rotation). The matter captured by the magnetic field of the neutron star at a given radius is orbiting with a given angular velocity, and will transfer angular momentum to the neutron star through the magnetic field lines. Outside $R_{co}$, this velocity is lower than the angular velocity of the neutron star, while it is higher inside. Therefore, roughly speaking, if $R_{in} < R_{co}$ the star spins up, and if $R_{in} > R_{co}$ it spins down. Various corrections can be made, integrating the torque from the matter outside and inside the corotation radius, and different authors have produced different prescriptions that can in general be treated by multiplying $R_{in}$ by a factor of order 1 (Ghosh & Lamb 1978; Wang 1996). When $R_{in} \sim R_{co}$, small changes of accretion rate move the inner radius back and forth around $R_{co}$, and we can expect the source to alternatively spin up and down. This situation is called spin equilibrium.

When $R_{in} \gg R_{co}$, the rotating magnetic field is able to swipe away the disk, and it is expected that accretion onto the neutron star will stop. This is known as propeller regime (Illarionov & Sunyaev 1975).

Around a super-Eddington accreting source, a third important radius is often cited, the spherization radius at which the disk departs from an ideal thin disk. Inside this radius, the mass in excess of the local Eddington limit is ejected in winds (Shakura & Sunyaev 1973):

$$R_{\text{sph}} = \frac{27}{4} R_g \frac{\dot{M}}{\dot{M}_{\text{Edd}}}, \qquad (C3)$$

where the Eddington mass-accretion rate $\dot{M}_{Edd} \approx L_{Edd}/\eta c^2 \approx 1.6 \cdot 10^{18} \text{ g s}^{-1}$ for a $1.4 M_{\odot}$ neutron star, where $\eta \approx 0.15$ is the efficiency, $R_g = GM/c^2$ is the gravitational radius, and c is the speed of light.

## Appendix D: Donor Star

The mass function determined through timing gives important insights on the kind of donor star we can expect:

$$f = \frac{M_d^3 \sin^3 i}{(M_p + M_d)^2} = \frac{\Omega_{\text{orb}}^2}{G} (a_p \sin i)^3 \approx 1.83 M_{\odot}, \qquad (D1)$$

where $\Omega_{orb} = 2\pi/P_{orb}$ is the orbital angular velocity, $a_p \sin i$ is the projected semimajor axis of the pulsar orbit, $M_d$ is the mass of the donor, $M_p$ is the mass of the pulsar, and $i$ is the inclination. In the formula above, the $\Omega_{orb}$ and $a_p \sin i$ are measured from pulsar timing, while the left-hand side can be used to infer the donor mass given reasonable assumptions about the pulsar mass and the inclination.

Since $\sin i$ cannot be larger than 1 (orbit edge-on), this poses a hard lower limit to the donor star mass that cannot be less than $3.56 M_{\odot}$ (assuming a neutron star mass of $1.4 M_{\odot}$). The absence of eclipses from a (most likely) Roche-lobe-filling donor pushes the lower limit to $\sim 5 M_{\odot}$ (B14) and corresponds to an upper limit on the inclination of $\sim 60^{\circ}$. An unlikely donor mass of $100 M_{\odot}$ corresponds instead to an inclination of $\sim 17^{\circ}$, which we take as a lower limit.

Similar arguments can be used to constrain the donor radius. Assuming Roche-lobe overflow, the size of the donor is fixed by the mass ratio and orbital separation.

With these constraints in mind (see Figure 4), and compared with known populations of donor stars in HMXBs, the most probable candidates are O/B giant stars between 5 and $100 M_{\odot}$. Between $\sim 17^{\circ}$ ($100 M_{\odot}$ donor) and $\sim 60^{\circ}$ ($5 M_{\odot}$ donor), we assume all orientations to be equally probable. This means that the values of the cosine of the inclination are equally probable between the two limiting cases $\cos 60^{\circ}$ and $\cos 17^{\circ}$. This gives an average inclination of $\sim 43^{\circ}$, corresponding to a donor mass of $\sim 7.8 M_{\odot}$. Note that an archival search in Hubble Space Telescope data found several stars of this range of masses which could in principle be the donor (Heida et al. 2019).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Roche-lobe radius in Eggleton approximation (Eggleton 1983) vs. mass for the donor star in M82 X-2. The donor has to lie around the black solid line in order to undergo Roche-lobe overflow. We overplot all donors from the HMXBs in Falanga et al. (2015) and all semidetached binary stars from Surkova & Svechnikov (2004) for comparison. For the HMXBs, we plot in black the Roche-lobe radius. Colors span the A (red) to O (blue) spectral types, and markers indicate different branches in the Hertzsprung-Russell diagram. The gray shaded area is excluded by the absence of eclipses. The black area is prohibited by the mass function and the necessity that $\sin i \le 1$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Neutron star dipolar magnetic field estimate assuming spin equilibrium as described in the text, comparing the models from Wang (1996; the area corresponds to values of $0.5 < \xi < 1$) and from Chashkina et al. (2017; the area covers the range of viscosity parameter $0.01 < \alpha_v < 0.3$), in four cases with a range of mass ratios and a different mass-loss fraction ($\alpha$ in Equation (B3)) from the donor. We show $R_{co}$ and $R_{sph}$ with vertical lines. Note that a lower mass ratio increases the estimated mass transfer. Also, the mass loss from the donor implies a larger mass transfer rate to justify the observed derivative. This explains the larger $R_{sph}$. We highlight the region where $0.8 < R_{in}/R_{co} < 1$, as expected from spin equilibrium. For the traditional model we took into account a change of accretion rate due to mass loss inside $R_{sph}$. The estimated dipolar magnetic field is always above $10^{13}$ G for reasonable values of the parameters.</figcaption></figure>

## Appendix E
### Magnetic Field Estimates

Traditional models, such as those proposed by Ghosh & Lamb (1978) or Wang (1996), consider a thin disk with negligible radiation effects, and the inner radius is given by Equation (C2). Therefore, given a mass-accretion rate, the position of the inner radius in this model is a function of the dipolar component of the magnetic field, modulo the order-unity constant $\xi$. Since the source is close to spin equilibrium, as demonstrated by the spin behavior over time (see Table 1), the inner radius has to be close to the corotation radius. Therefore, equating the inner radius to the corotation radius, we can get an estimate of the magnetic field, as shown in Figure 5 with the blue band. Here, we take into account mass losses in a wind inside $R_{sph}$, when relevant, until $R_{in}$. We use the “classical” mass loss obtained when the effects of advection are neglected (Shakura & Sunyaev 1973). In this case, the accretion rate drops linearly with radius, and thus an upper limit on the accretion rate (i.e., an upper limit on the mass-loss rate) corresponds to an upper limit on $R_{in}$ and a lower limit to the magnetic field strength. Despite this conservative approach, the estimate on the magnetic field is robustly above $10^{13}$ G. Most models for sub-Eddington accretion agree within an order of magnitude for the treatment of spin equilibrium (see, e.g., Chen et al. 2021). However, it is possible that these models, which are based on the interaction of a magnetized neutron star with a thin disk, with no radiation pressure either from the disk or from the central object, need to be corrected in the case of

super-Eddington disks. Chashkina et al. (2017, 2019) have investigated this issue, finding that, indeed, the disk structure changes significantly when radiation pressure becomes dominant. In particular, they find that $\xi$ is not constant, but depends on local (inside the disk) and external (e.g., from the neutron star) radiation pressure, and the amount of advection in the disk. With the transfer rate $>100 \times$ Eddington we infer in this paper, the inner radius becomes almost independent of the mass-accretion rate and is described by Equation (61) from Chashkina et al. (2017):

$$ \frac{R_{in}}{R_g} \approx \left( \frac{73\alpha}{24} \right)^{2/9} \left[ \lambda \left( \frac{\mu}{10^{30} \text{G cm}^3} \right)^2 \right]^{2/9}, \quad (E1) $$

where $\alpha \sim 0.1$ is the viscosity in the disk and $\lambda \sim 4 \cdot 10^{10} (M_p/1.4M_\odot)^{-5}$. Figure 5 shows that, for a reasonable range of the viscosity parameter<fnref n="19" /> $0.01 < \alpha_v < 0.3$, the estimate of the magnetic field obtained by equating the inner radius to the corotation radius using Equation (E1) is similar to the prediction of traditional models using Equation (C2), confirming an estimated magnetic field for M82 X-2 above $10^{13}$ G, as estimated with the classical model and by other authors in the literature (Tsygankov et al. 2016; Chen et al. 2021).

<footnote n="19">¹⁹ We call this parameter $\alpha_v$ instead of $\alpha$ to avoid confusion with the mass-loss parameter $\alpha$.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Pulsed amplitude vs. energy, using XMM-Newton data from ObsIDs 0870940101 and 0870940401, and NuSTAR data from ObsIDs 30602027002 and 30602027004. We overplot the raw count spectra of NuSTAR and XMM, showing how XMM-Newton has many more counts, but in energy intervals with little or no pulsations.</figcaption></figure>

## Appendix F
### Pulsed Fraction in the XMM-Newton Energy Band

As opposed to many other PULXs, M82 X-2 is very difficult to study with XMM-Newton. Pulsations were not detected in many past observations of M82 X-2, despite the higher angular resolution of the EPIC-pn instrument. One of the reasons is the lack of pulsations below 3 keV, due to both an intrinsic low pulsed amplitude and the very strong emission of M82 X-1 and the M82 galaxy itself that increase the background at low energies. In the 2021 quasi-simultaneous observations with XMM-Newton and NuSTAR, we did manage to detect pulsations with EPIC-pn (Figure 6). M82 X-2 was observed by XMM-Newton on UT 2021-04-06 and 2021-04-16 for a total on-source exposure time of ~70 ks. The only camera on board XMM-Newton that is able to detect pulsations from M82 X-2 is EPIC-pn, which was set in Full Window mode.

We downloaded the data from the two observations from the XMM-Newton archive<fnref n="20" /> and processed them with the Science Analysis Software (SAS) version 20211130.

We ran the standard pipeline, using the tool *epchain* to obtain cleaned event files. The M82 field is very crowded, and it is not possible to separate the emission of M82 X-2, M82 X-1, and the diffuse galactic center emission. However, being mostly interested in the timing properties of M82 X-2, a precise modeling of the background is not strictly needed. We selected photons coming from a region of ~50" around the putative position of M82 X-2. We cleaned the data from periods of high background activity. Finally, we barycentered the data using the tool *barycen* using the Chandra position of M82 X-2, with the same ephemeris used in *barycorr*.

After this preprocessing, we folded the cleaned and barycentered event lists at the ephemeris obtained from the nearest NuSTAR observations, slightly adjusting the spin frequency through the maximization of the Rayleigh test. We calculated the pulsed fraction from a sinusoidal modeling of the pulsed profile, as $(Max-Min)/(Max + Min)$. We plot this

pulsed fraction, and the corresponding pulsed fraction from the quasi-simultaneous NuSTAR observations, in Figure 6.

## ORCID iDs

Matteo Bachetti https://orcid.org/0000-0002-4576-9337

Marianne Heida https://orcid.org/0000-0002-1082-7496

Thomas Maccarone https://orcid.org/0000-0003-0976-4755

Daniela Huppenkothen https://orcid.org/0000-0002-1169-7486

Gian Luca Israel https://orcid.org/0000-0001-5480-6438

Didier Barret https://orcid.org/0000-0002-0393-9190

Murray Brightman https://orcid.org/0000-0002-8147-2602

McKinley Brumback https://orcid.org/0000-0002-4024-6967

Hannah P. Earnshaw https://orcid.org/0000-0001-5857-5622

Karl Forster https://orcid.org/0000-0001-5800-5531

Felix Fürst https://orcid.org/0000-0003-0388-0560

Brian W. Grefenstette https://orcid.org/0000-0002-1984-2932

Amruta D. Jaodand https://orcid.org/0000-0002-3850-6651

Kristin K. Madsen https://orcid.org/0000-0003-1252-4891

Matthew Middleton https://orcid.org/0000-0002-8183-2970

Sean N. Pike https://orcid.org/0000-0002-8403-0041

Maura Pilia https://orcid.org/0000-0001-7397-8091

Juri Poutanen https://orcid.org/0000-0002-0983-0049

Daniel Stern https://orcid.org/0000-0003-2686-9241

John A. Tomsick https://orcid.org/0000-0001-5506-9855

Dominic J. Walton https://orcid.org/0000-0001-5819-3552

Jörn Wilms https://orcid.org/0000-0003-2065-5410

## References

Abarca, D., Parfrey, K., & Kluźniak, W. 2021, ApJ, 917, L31

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Bachetti, M. 2018, HENDRICS: High Energy Data Reduction Interface from the Command Shell, Astrophysics Source Code Library, ascl:1805.019

Bachetti, M., Pilia, M., Huppenkothen, D., et al. 2021, ApJ, 909, 33

Bachetti, M., Rana, V., Walton, D. J., et al. 2013, ApJ, 778, 163

<footnote n="20">²⁰ nxsa.esac.esa.int</footnote>
Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nat, 514, 202

Bachetti, M., Maccarone, T. J., Brightman, M., et al. 2020, ApJ, 891, 44

Barcons, X., Barret, D., Decourchelle, A., et al. 2017, AN, 338, 153

Blundell, K. M., Bowler, M. G., & Schmidtobreick, L. 2008, ApJL, 678, L47

Brice, N., Zane, S., Turolla, R., & Wu, K. 2021, MNRAS, 504, 701

Brightman, M., Harrison, F. A., Bachetti, M., et al. 2019, ApJ, 873, 115

Buccheri, R., Bennett, K., Bignami, G. F., et al. 1983, A&A, 128, 245

Chashkina, A., Abolmasov, P., & Poutanen, J. 2017, MNRAS, 470, 2799

Chashkina, A., Lipunova, G., Abolmasov, P., & Poutanen, J. 2019, A&A, 626, A18

Chen, X., Wang, W., & Tong, H. 2021, JHEAp, 31, 1

Cherepashchuk, A. M., Belinski, A. A., Dodin, A. V., & Postnov, K. A. 2021, MNRAS, 507, L19

Chernov, S. V. 2020, ARep, 64, 425

Dall'Osso, S., Perna, R., & Stella, L. 2015, MNRAS, 449, 2144

D'Angelo, C. R., & Spruit, H. C. 2012, MNRAS, 420, 416

Eggleton, P. P. 1983, ApJ, 268, 368

Ekşi, K. Y., Andaç, C., Çıkintoğlu, S., et al. 2015, MNRAS Let, 448, L40

Erkut, M. H., Türkoğlu, M. M., Ekşi, K. Y., & Alpar, M. A. 2020, ApJ, 899, 97

Fabrika, S., Karpov, S., Abolmasov, P., Sholukhova, O., & Fabbiano, G. 2006, in Populations of High Energy Sources in Galaxies Proc. of the 230th Symp. of the IAU (Cambridge: Cambridge Univ. Press), 278

Fabrika, S. N., Atapin, K. E., Vinokurov, A. S., & Sholukhova, O. N. 2021, AstBu, 76, 6

Falanga, M., Bozzo, E., Lutovinov, A., et al. 2015, A&A, 577, A130

Foreman-Mackey, D. 2016, JOSS, 1, 24

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition (Accretion Power in Astrophysics) (Cambridge: Cambridge Univ. Press)

Ghosh, P., & Lamb, F. K. 1978, ApJL, 223, L83

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

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Heida, M., Harrison, F. A., Brightman, M., et al. 2019, ApJ, 871, 231

Hunter, J. D. 2007, CSE, 9, 90

Huppenkothen, D., Bachetti, M., Stevens, A. L., Migliari, S., & Balm, P. 2016, Stingray: Spectral-timing software, Astrophysics Source Code Library, ascl:1608.001

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

Illarionov, A. F., & Sunyaev, R. A. 1975, A&A, 39, 185

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

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

Joss, P. C., & Rappaport, S. A. 1984, ARA&A, 22, 537

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

Kelley, R., Rappaport, S., & Petre, R. 1980, ApJ, 238, 699

King, A., & Lasota, J.-P. 2020, MNRAS, 494, 3611

King, A., Lasota, J.-P., & Kluźniak, W. 2017, MNRAS Let., 468, L59

King, A. R. 2008, MNRAS Let., 385, L113

Kosec, P., Pinto, C., Walton, D. J., et al. 2018, MNRAS, 479, 3978

Lam, S. K., Pitrou, A., & Seibert, S. 2015, Proc. of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM '15 (New York: ACM), doi:10.1145/2833157.2833162

Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274

Levine, A., Rappaport, S., Deeter, J. E., Boynton, P. E., & Nagase, F. 1993, ApJ, 410, 328

Levine, A. M., Rappaport, S. A., & Zojcheski, G. 2000, ApJ, 541, 194

Lu, W., Fuller, J., Quataert, E., & Bonnerot, C. 2022, arXiv:2204.00847

Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

Madsen, K., Hickox, R., Bachetti, M., et al. 2019, BAAS, 51, 166

Mellah, I. E., Sundqvist, J. O., & Keppens, R. 2019, A&A, 622, L3

Middleton, M. J., Higginbottom, N., Knigge, C., Khan, N., & Wiktorowicz, G. 2022, MNRAS, 509, 1119

Middleton, M. J., Walton, D. J., Alston, W., et al. 2021, MNRAS, 506, 1045

Misra, D., Fragos, T., Tauris, T. M., Zapartas, E., & Aguilera-Dena, D. R. 2020, A&A, 642, A174

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

Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015, MNRAS, 454, 2539

Pinto, C., Fabian, A., Middleton, M., & Walton, D. 2016, AN, 338, 234

Pletsch, H. J., & Clark, C. J. 2015, ApJ, 807, 18

Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187

Ray, P. S., Arzoumanian, Z., Brandt, S., et al. 2019, BAAS, 51, 231

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

Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2004, ApJL, 616, L151

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

Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327, 620

Stoyanov, K. A., & Zamanov, R. K. 2009, AN, 330, 727

Surkova, L. P., & Svechnikov, M. A. 2004, VizieR On-line Data Catalog: V/115

Tauris, T. M., & Savonije, G. J. 2001, in The neutron star - black hole connection: Proc. of the NATO Advanced Study Institute, ed. C. Kouveliotou, J. Ventura, & E. Van den Heuvel, Vol. 567 (Dordrecht: Kluwer), 337

Tauris, T. M., & van den Heuvel, E. P. J. 2006, in Compact Stellar X-Ray Sources, ed. W. Lewin & M. van der Klis (Cambridge: Cambridge Univ. Press), 623

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

van den Heuvel, E. P. J. 1994, Saas-Fee Advanced Course of the Swiss Society for Astrophysics and Astronomy: Interacting binaries (Berlin: Springer), 263

Vasilopoulos, G., Koliopanos, F., Haberl, F., et al. 2021, ApJ, 909, 50

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Walton, D. J., Tomsick, J. A., Madsen, K. K., et al. 2016, ApJ, 826, 87

Wang, Y.-M. 1996, ApJL, 465, L111

Xu, K., & Li, X.-D. 2017, ApJ, 838, 98

Zhang, S. N., Feroci, M., Santangelo, A., et al. 2017, Proc. SPIE, 9905, 99051Q