# All at Once: Transient Pulsations, Spin-down, and a Glitch from the Pulsating Ultraluminous X-Ray Source M82 X-2

Matteo Bachetti<sup>1,2,12</sup> (✉), Thomas J. Maccarone<sup>3</sup>, Murray Brightman<sup>2</sup>, McKinley C. Brumback<sup>4</sup> (✉), Felix Fürst<sup>5</sup> (✉), Fiona A. Harrison<sup>2</sup>, Marianne Heida<sup>2,6</sup> (✉), Gian Luca Israel<sup>7</sup> (✉), Matthew J. Middleton<sup>8</sup> (✉), John A. Tomsick<sup>9</sup> (✉), Natalie A. Webb<sup>10</sup>, and Dominic J. Walton<sup>11</sup> (✉)

## Abstract

The first pulsating ultraluminous X-ray source (PULX) to be identified is M82 X-2. After the discovery in 2014, NuSTAR observed the M82 field 15 times throughout 2015 and 2016. In this paper, we report the results of pulsation searches in all of these data sets and find only one new detection. This new detection allows us to refine the orbital period of the source and measure an average spin-down rate between 2014 and 2016 of ~−6 × 10<sup>−11</sup> Hz s<sup>−1</sup>, which is in contrast to the strong spin-up seen during the 2014 observations, representing the first detection of spin-down in a PULX system. Thanks to the improved orbital solution allowed by this new detection, we are also able to detect pulsations in additional segments of the original 2014 data set. We find a glitch superimposed on the very strong and variable spin-up already reported—the first positive glitch identified in a PULX system. We discuss the new findings in the context of current leading models for PULXs.

*Unified Astronomy Thesaurus concepts: Pulsars* (1306); *Pulsar timing method* (1305); *Accretion* (14); *X-ray astronomy* (1810); *Orbital motion* (1179); *Time series analysis* (2109); *Fast Fourier transform* (1958); *Neutron stars* (1108); *Compact objects* (288)

## 1. Introduction

Ultraluminous X-ray sources (ULXs) are off-nuclear point sources with X-ray luminosities exceeding the Eddington limit for a stellar-remnant black hole (see Kaaret et al. 2017, for a review). In 2014, Bachetti et al. (2014, hereafter B14) reported pulsations from the known ULX M82 X-2 showing that at least some ULXs are neutron stars. With additional observations and timing analysis, more of these objects are being found to pulsate (e.g., NGC 5907 X-1: Israel et al. 2017a; NGC 7793 P13: Israel et al. 2017b, Fürst et al. 2016; NGC 300 X-1: Carpano et al. 2018; NGC 1313 X-2: Sathyaprakash et al. 2019). The spectral and variability properties of these pulsating ULXs (PULXs<sup>13</sup>) are similar to the bulk of the ULX population. Looking in detail, they tend to have slightly harder spectra and a higher level of long-term variability than average.<sup>14</sup> This fact, together with the intrinsic difficulties in finding pulsations in these distant sources, means that it is plausible that most of the ULX population is powered by neutron stars.<sup>15</sup>

The PULXs found so far have in common periods around 1 s (with the notable exception of NGC 300 X-1, which is much

slower) and strong spin-up during and between observations spaced by months to years. The spin-up is likely due to the torque of an accretion disk onto the star as observed in many accreting neutron stars in X-ray binary systems (e.g., Bildsten & Brown 1997; Rappaport & Joss 1977; Burderi et al. 2006; Perna et al. 2006).

Accretion onto magnetized neutron stars happens through the interaction of the accretion disk with the strong magnetic field. The general framework of accretion onto magnetized objects is described in a number of papers (e.g., Ghosh & Lamb 1979a, 1979b; Wang 1995) and reviewed in Frank et al. (2002). The disk is truncated at a distance $R_{in}$ (often called the Alfvén radius) from the neutron star, where magnetic stresses equal the ram pressure from the matter in the disk. Therefore, the same $R_{in}$ can be determined by hugely different magnetic field strengths, provided that they are balanced by different mass accretion rates. The magnetic field, besides interrupting the disk at $R_{in}$, also penetrates the disk beyond the inner radius. The interaction between the matter and the magnetic field at a given distance from the neutron star produces a torque on the star that depends on the relative angular velocity of the matter with respect to the star. Calling $R_{co}$ the corotation radius, the radius at which the Keplerian angular velocity of the accreting matter equals the angular velocity of the star, the total torque depends on how many field lines thread the disk inside and beyond $R_{co}$ (Wang 1995). The observation of spin-up and spin-down in accreting pulsars with plausibly very different magnetic fields and companion stars is usually explained

<footnote n="1"><sup>1</sup> INAF-Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius, Italy; matteo.bachetti@inaf.it</footnote>

<footnote n="2"><sup>2</sup> Space Radiation Laboratory, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA</footnote>

<footnote n="3"><sup>3</sup> Department of Physics and Astronomy, Texas Tech University, Lubbock, TX, USA</footnote>

<footnote n="4"><sup>4</sup> Department of Physics & Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH 03755, USA</footnote>

<footnote n="5"><sup>5</sup> European Space Astronomy Centre (ESA/ESAC), Operations Department, Villanueva de la Canada Madrid, Spain</footnote>

<footnote n="6"><sup>6</sup>European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="7"><sup>7</sup> Osservatorio Astronomico di Roma, INAF, via Frascati 33, I-00078 Monte Porzio Catone, Italy</footnote>

<footnote n="8"><sup>8</sup> Department of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, UK</footnote>

<footnote n="9"><sup>9</sup> Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA</footnote>

<footnote n="10"><sup>10</sup> IRAP, Université de Toulouse, CNRS, UPS, CNES, 9 Avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="11"><sup>11</sup> Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote>Received 2019 May 15; revised 2020 January 11; accepted 2020 January 15; published 2020 March 3</footnote>

<footnote n="12"><sup>12</sup> Fulbright Visiting Scholar.</footnote>

<footnote n="13"><sup>13</sup> Also referred to as ultraluminous X-ray pulsars.</footnote>

<footnote n="14"><sup>14</sup> For a detailed view, see Koliopanos et al. (2017), Pintore et al. (2017) and Walton et al. (2018).</footnote>

<footnote n="15"><sup>15</sup> For deeper investigations of ULX populations, see Middleton & King (2017) and Wiktorowicz et al. (2018).</footnote>
in this theoretical framework with good success. As a final ingredient in standard accreting pulsars, the X-ray luminosity is considered a proxy of the mass accretion rate (following Shakura & Sunyaev 1973), as for thin disks, we expect a relation $L \propto \dot{M}$.

However, the luminosity of PULXs is far higher than the Eddington luminosity for a neutron star, even considering the modifications to the Eddington limit that arise in pulsar accretion columns (e.g., Basko & Sunyaev 1976). This has led to a lively theoretical debate, which is still ongoing.

Because of this super-Eddington accretion, any treatment of PULXs should include a third special radius, called the spherization radius, $R_{\text{sph}}$. In super-Eddington disks, inside $R_{\text{sph}}$, we cannot assume a standard thin disk configuration. Inside this radius, the disk is swollen and strong winds are ejected due to radiation pressure, limiting the total amount of matter that actually approaches the inner radius and participates in the acceleration of the star (Shakura & Sunyaev 1973; Poutanen et al. 2007). Moreover, the observed luminosity itself can be misleading. The bolometric luminosity of a super-Eddington disk increases according to

$$L_{\text{acc}} \approx L_{\text{Edd}}(1 + \ln \dot{m}), \quad (1)$$

where $\dot{m} = \dot{M}/\dot{M}_{\text{Edd}}$. However, this is probably not the bolometric luminosity measured by the observer, because radiation might be beamed. King (2009) described a geometry where the origin of the beaming is the wind. In their model, the wind angle is proportional to the mass accretion rate, and this leads to a beaming factor $b \propto \dot{M}^{-2}$. In the end, one obtains

$$L_{\text{obs}} = L_{\text{acc}}/b \propto L_{\text{Edd}}\dot{m}^2(1 + \ln \dot{m}). \quad (2)$$

With this in mind, the literature about PULXs has explored two main cases: $R_{\text{sph}} \gtrsim R_{\text{in}}$ and $R_{\text{sph}} < R_{\text{in}}$ (for recent examples, see Walton et al. 2018; Middleton et al. 2019a; Mushtukov et al. 2019).

If $R_{\text{sph}} > R_{\text{in}}$, the disk becomes locally super-Eddington before reaching the inner radius. So, the actual mass accretion rate (the mass that reaches $R_{\text{in}}$ and accretes onto the neutron star) is lower than inferred from the luminosity because winds eject a large fraction of the matter, and, additionally, the luminosity itself is beamed and overestimated. The observation of blueshifted lines reminiscent of strong winds and various kinds of nebulae around ULXs supports this hypothesis (e.g., Grisé et al. 2011; Pinto et al. 2017). However, the bulk of observations of nebulae, wind signatures, and broad pulse profiles exclude extreme beaming factors (see review by Kaaret et al. 2017).

In the opposite case, the disk remains thin until it reaches $R_{\text{in}}$, and no strong winds are launched. The effects of the huge accretion rate manifest themselves in the accretion column, giving rise to an optically thick envelope around the star that can, in principle, suppress high-frequency variability (Mushtukov et al. 2017). Since the inferred $R_{\text{in}}$ is quite large and the inferred magnetic fields are magnetar-like, this model justifies the high luminosity because, in very high magnetic fields, the Thompson scattering cross section—one of the key ingredients of the Eddington limit—is modified, and the Eddington limit itself is much higher (Basko & Sunyaev 1976; Mushtukov et al. 2015). In this interpretation, all discovered PULXs have rather slow pulsations because pulsations with higher frequencies are washed out by the envelope.

One way to address this controversy is to try to determine the mass accretion rate in an independent way and use it to calculate the spherization radius. Since PULXs almost certainly

have super-Eddington mass accretion rates, whether or not this material ends up on the compact object, we expect that this large mass transfer produces visible effects on the binary system. Motivated by this, we obtained a long NuSTAR observation of the pulsar with the aim of detecting pulsations and, through the measurement of a possible orbital shrinking, constraining the total mass exchange in the system. The observation was performed on UT 2016 April 15–19 (MJD 57,493.29–57,497.34).

NuSTAR does not resolve M82 X-2 from another ULX, M82 X-1, only 5" away. This source is known to reach luminosities an order of magnitude higher than M82 X-2, and only Chandra is able to resolve the two ULXs. Therefore, we undertook another program to monitor M82 monthly with Chandra with 25 ks per pointing with simultaneous NuSTAR 40 ks pointings. This Chandra–NuSTAR program was designed to study the spectral evolution of the two sources, as well as to characterize the overall binary population in the M82 galaxy (see Brightman et al. 2019, 2020).

In this paper, we present the timing analysis of the new data listed above and show that the pulsations are only detected in one new observation, despite the fact that, in some cases, we infer a high enough flux level from M82 X-2 relative to M82 X-1 from the Chandra observations for the pulsations to be detected. Revisiting the work by B14, we show that the pulsed fraction evolves independently from the flux from the nearby M82 X-1, showing that this behavior is intrinsic to M82 X-2. Thanks to the new detection, we are able to measure the orbital period of the system with greater precision, recover pulsations in one more old observation, and detect a pulsar glitch.

In Section 2 we describe the data reduction procedure. In Section 3 we detail the pulsation searches performed and report on the new pulsations and the orbit and spin measurements. We discuss the results in Section 4 and summarize them in Section 5.

## 2. Data Reduction

The observations considered in this work are detailed in Table 1. For NuSTAR, we used data processed with nupipeline from HEASOFT v.6.25 (HEASARC 2014), with standard options. We barycentered old and new data in two independent ways. First, we used the FTOOL barycorr (Blackburn 1995) with the standard CALDB clock file, then we tested a new temperature-driven model of the clock offsets (under development)<fnref n="16" /> and used PINT<fnref n="17" /> for the satellite orbit calculations. The pulse period of M82 X-2 is 3 orders of magnitude longer than the observed difference between the two barycentering methods. For the orbit file, we used the attitude-orbit file produced by nupipeline.<fnref n="18" /> Running the pipeline using default values left some intervals of increased background activity (probably due to the South Atlantic Anomaly), producing spikes in the light curve. We verified that that the spikes correspond to increased activity in the entire field of view, not just from a single source; we removed the relevant intervals by modifying the good time intervals (GTIs). Following B14, we selected photons from a

<footnote n="16">¹⁶ The temperature-driven clock correction is now used to produce the official clock correction files with NuSTAR. The new clock correction is the default since barycorr v.2.2. M. Bachetti & C. B. Markwardt et al. 2020, in preparation.</footnote>

<footnote n="17">¹⁷ www.github.com/nanograv/PINT</footnote>

<footnote n="18">¹⁸ During the study of the temperature-driven model for the spacecraft clock, we realized that the orbit file distributed in the auxiliary data of NuSTAR observations and recommended for use with barycorr did not properly account for leap seconds.</footnote>
<table><caption>**Table 1**
Pulsation Search Results from All Available Observations of M82 X-2</caption><thead><tr><th>ObsID</th><th>MJD</th><th>Date</th><th>Length<sup>a</sup><br>(ks)</th><th>Exposure<sup>a</sup><br>(ks)</th><th>Count Rate<sup>b</sup><br>(s<sup>-1</sup>)</th><th>Meas. p.f.<sup>c</sup><br>ΔF<sub>X2</sub>/F<sub>tot</sub>(%)</th><th>Flux Ratio<sup>d</sup><br>F<sub>X2</sub>/F<sub>X1</sub>(%)</th><th>Corr. p.f.<sup>c</sup><br>ΔF<sub>X2</sub>/F<sub>X2</sub>(%)</th><th>L<sub>39</sub></th></tr></thead><tbody><tr><td>80002092002</td><td>56,681</td><td>2014 Jan 23</td><td>123</td><td>66</td><td>0.9</td><td>&lt;5</td><td>...</td><td>...</td><td>...</td></tr><tr><td>80002092004<sup>f</sup></td><td>56,683</td><td>2014 Jan 25</td><td>171</td><td>90</td><td>1.0</td><td>4</td><td>...</td><td>...</td><td>...</td></tr><tr><td>80002092006</td><td>56,686</td><td>2014 Jan 28</td><td>579</td><td>310</td><td>0.9</td><td>5</td><td>...</td><td>...</td><td>...</td></tr><tr><td>80002092007</td><td>56,692</td><td>2014 Feb 4</td><td>562</td><td>306</td><td>0.9</td><td>7</td><td>132*</td><td>12</td><td>8.1*</td></tr><tr><td>80002092008</td><td>56,699</td><td>2014 Feb 10</td><td>62</td><td>34</td><td>1.0</td><td>7</td><td>...</td><td>...</td><td>...</td></tr><tr><td>80002092009</td><td>56,700</td><td>2014 Feb 11</td><td>213</td><td>115</td><td>0.9</td><td>9</td><td>...</td><td>...</td><td>...</td></tr><tr><td>80002092011</td><td>56,720</td><td>2014 Mar 3</td><td>201</td><td>111</td><td>0.6</td><td>3</td><td>...</td><td>...</td><td>...</td></tr><tr><td>50002019002</td><td>57,038</td><td>2015 Jan 15</td><td>56</td><td>31</td><td>0.7</td><td>&lt;6</td><td>40</td><td>&lt;27</td><td>5.7</td></tr><tr><td>50002019004</td><td>57,042</td><td>2015 Jan 19</td><td>283</td><td>161</td><td>0.7</td><td>&lt;4</td><td>33*</td><td>&lt;14</td><td>2.0*</td></tr><tr><td>90101005002</td><td>57,194</td><td>2015 Jun 20</td><td>56</td><td>37</td><td>2.3</td><td>&lt;3</td><td>32</td><td>&lt;17</td><td>19.8</td></tr><tr><td>80202020002</td><td>57,414</td><td>2016 Jan 26</td><td>66</td><td>36</td><td>1.5</td><td>&lt;4</td><td>0*</td><td>N.A.</td><td>0.0*</td></tr><tr><td>80202020004</td><td>57,442</td><td>2016 Feb 23</td><td>61</td><td>32</td><td>1.3</td><td>&lt;5</td><td>52</td><td>&lt;18</td><td>13.9</td></tr><tr><td>80202020006</td><td>57,483</td><td>2016 Apr 5</td><td>54</td><td>31</td><td>0.8</td><td>&lt;6</td><td>44</td><td>&lt;26</td><td>6.7</td></tr><tr><td>30101045002</td><td>57,493</td><td>2016 Apr 15</td><td>350</td><td>189</td><td>1.0</td><td>&lt;3</td><td>...</td><td>...</td><td>...</td></tr><tr><td>80202020008</td><td>57,503</td><td>2016 Apr 24</td><td>67</td><td>40</td><td>1.1</td><td>&lt;5</td><td>79</td><td>&lt;13</td><td>10.8</td></tr><tr><td>30202022002</td><td>57,543</td><td>2016 Jun 3</td><td>60</td><td>39</td><td>1.0</td><td>&lt;5</td><td>0*</td><td>N.A.</td><td>0.1*</td></tr><tr><td>30202022004</td><td>57,571</td><td>2016 Jul 1</td><td>68</td><td>47</td><td>1.7</td><td>&lt;4</td><td>33</td><td>&lt;17</td><td>15.1</td></tr><tr><td>30202022008</td><td>57,599</td><td>2016 Jul 29</td><td>67</td><td>43</td><td>1.6</td><td>&lt;4</td><td>4*</td><td>N.A.</td><td>1.2*</td></tr><tr><td>30202022010</td><td>57,619</td><td>2016 Aug 19</td><td>69</td><td>44</td><td>1.2</td><td>&lt;4</td><td>16*</td><td>&lt;37</td><td>4.3*</td></tr><tr><td>90201037002</td><td>57,641</td><td>2016 Sep 10</td><td>94</td><td>80</td><td>1.2</td><td>3</td><td>...</td><td>...</td><td>...</td></tr><tr><td>90202038002</td><td>57,669</td><td>2016 Oct 7</td><td>71</td><td>45</td><td>0.9</td><td>&lt;5</td><td>5*</td><td>N.A.</td><td>0.8*</td></tr><tr><td>90202038004</td><td>57,723</td><td>2016 Nov 30</td><td>68</td><td>43</td><td>0.8</td><td>&lt;5</td><td>0*</td><td>N.A.</td><td>0.0*</td></tr></tbody><tfoot><tr><td colspan="10">Det. puls.</td><td>1000 (58%)</td></tr><tr><td colspan="10">Undet.</td><td>730 (42%)</td></tr><tr><td colspan="10">Low flux</td><td>207</td></tr></tfoot></table>

**Notes.** Intervals where M82 X-2 is weaker than 10% M82 X-1 are not counted in the statistics of nondetections.

<sup>a</sup> The length of the observation is the UT stop time minus the UT start time. The exposure is the actual on-source time not including occultation from the Earth and other “bad” intervals.

<sup>b</sup> Count rates using FPMA + FPMB.

<sup>c</sup> The pulsed fraction is referred to the total X-ray flux of M82, since M82 X-1 and X-2 are not separable in *NuSTAR*.

<sup>d</sup> The flux ratio between M82 X-2 and M82 X-1 is estimated from *Chandra* data, either from the detailed and low-pileup spectral modeling of Brightman et al. (2020) or from Brightman et al. (2019; indicated by an asterisk). Ellipses indicate that a *Chandra* data set was not available.

<sup>e</sup> A corrected estimate of the pulsed fraction when a flux ratio is available. Not allowed (N.A.) indicates that the pulsed flux should be higher than the flux measured by *Chandra*.

<sup>f</sup> Detected only after MJD 56,683.5.

region of 70" around M82 X-2. This region contains a large number of X-ray sources (e.g., B14; Brightman et al. 2016), with the total flux typically dominated by the combination of M82 X-1 ($L_{x,max} \sim 10^{41}$ erg s$^{-1}$) and M82 X-2 ($L_{x,max} \sim 2 \times 10^{40}$ erg s$^{-1}$; Brightman et al. 2019). Given their $\sim 5''$ separation, it is not possible to disentangle the contribution from these two sources to the *NuSTAR* data, so they must be considered together.

We did not redo the *Chandra* data reduction, and we used the data from Brightman et al. (2019, 2020).

### 3. Search for Pulsations and Results

We used a number of methods to detect new pulsations from M82 X-2, ranging from a focused search around the known frequency values and orbital parameters to a quasi-blind search (allowing for very large variations of the pulse frequency).

Data analysis was done through a set of custom Python scripts based on HENDRICS<fnref n="19" /> (Bachetti 2018), Stingray<fnref n="20" />

(Huppenkothen et al. 2016, 2019), PINT (Luo et al. 2019), and PRESTO (Ransom 2001, 2011).

In the 2 Ms long campaign performed in 2014, where the pulsar was first detected, the pulsations did not have a constant rms. A large variation of the pulsed fraction was reported by B14 and is also shown in Figure 4. Orbital motion represented an additional difficulty. In fact, in 2014, the pulsar was initially detected only after the pulsed fraction had increased well above the detection level, because orbital motion smears out the observed pulsed frequency, and the fainter signal was dominated by white noise. Only after the first detection was made at a high pulsed fraction and the orbital parameters were determined was the pulse found in some of the previous observations.

Therefore, any search for pulsations in new observations had to account for at least three complications: weak pulsations, strong spin-up, and orbital motion. We used a two-tiered approach to this search: a deep search for pulsations around the expected spin period values and orbital parameters and a more general search using multiple spin derivatives on 40 ks segments of data.

We report the details in the following subsections.

<footnote n="19">¹⁹ www.github.com/StingraySoftware/HENDRICS</footnote>

<footnote n="20">²⁰ www.github.com/StingraySoftware/stingray</footnote>
## 3.1. Deep Pulsation Searches—First Pass

The first attempt consisted of a deep search of pulsations using the $Z_2^2$ statistics (Buccheri et al. 1983). The $Z_2^2$ statistics was calculated from the functions in stingray and HENDRICS.<fnref n="21" />

We varied the frequency between the observed frequency in 2014 and the maximum frequency expected from a constant source spin-up of $5 \times 10^{-10}$ Hz s$^{-1}$ (more than double the maximum spin-up observed in 2014). The choice of the frequency interval was driven primarily by computing time and using an acceptable number of trial values. The frequency step was eight times finer than the standard (e.g., fast Fourier transform (FFT)) $\delta\nu = 1/T$, with T the observing time. This was to avoid the effects of spectral leakage on weak pulsations (e.g., see discussion of “interbinning” by Ransom et al. 2002). Not finding new detections above the standard $3 - \sigma$ detection level accounting for the number of trials, we shifted the orbital phase by trial values spaced by ~200 s in order to account for a possibly imperfect orbital solution.<fnref n="22" /> Again (and taking into account the increased number of trials), we did not find new significant pulsations with this strategy.

## 3.2. Accelerated Search

We used the PRESTO suite of pulsar search programs (Ransom 2001; Ransom et al. 2002) to run two different techniques of pulsation search, one Fourier-based (with the tool accelsearch) and one epoch folding-based (with the tool prepfold).

Following the strategy adopted in 2014, we split the light curves into segments of 30 ks, with sliding windows overlapping by a factor of ~0.5. This is motivated by the following: a longer light curve, in principle, allows for more signal to be accumulated in the periodograms, but orbital motion smears the signal of an accelerated search if the length is more than a certain fraction of the orbital period (the usual rule of thumb is $1/7$ of $P_{\text{orb}}$). This was clearly seen in the 2014 campaign, where the maximum detectability with the accelerated methods in prepfold was indeed obtained with segments of 30 ks, or about $1/7P_{\text{orb}}$.

We binned the light curves to 1 ms and produced binary floating point data sets in a format understandable by PRESTO using the HENbinary script in HENDRICS.

We first used the accelsearch tool, which performs an accelerated search of pulsations based on the FFT and matched filters. We used the lowest detection limit (-$\sigma_{\text{lim}}=1$ on the command line) and specified that the data were obtained by photons (-photon). A number of “candidates” (possible signals above the threshold power value) were produced by this search and followed up with prepfold to evaluate the significance using epoch folding. All candidate periods from this search were very different from the reasonable interval of pulsations from M82 X-2. We could also safely dismiss the possibility that these candidates were from other pulsars in the field of view, as the candidate pulsations always had very low significance and were not consistent between segments of data,

as one would expect from pulsar candidates. Some of these candidates are suspiciously close to the beats of the fundamental frequencies of the NuSTAR detector<fnref n="23" /> or the orbital data gaps.

We next tried to use the tool prepfold to run an epoch folding-based search (i.e., suggesting candidate periods instead of letting accelsearch do so). The prepfold tool automatically searches the $p-\dot{p}$ plane, including all reasonable $\dot{p}$ values expected from the source if the interval length is 30 ks. Even if the maximum $\dot{p}$ in 2014 was high, many accreting pulsars have an average $\dot{p}$ that is much smaller than the instantaneous value found in single observations due to the alternation of spin-up and spin-down events that the accretion torque produces when the accretion rate increases and decreases. Therefore, we specified a starting value of the period around the mean value in 2014 and then used a number of starting values further and further from the mean value in both directions, randomly distributed in an interval larger than 10 times the period variation expected if the maximum $\dot{p}$ in 2014 was constant until our new observations (which would result in spin periods of 1.31–1.372 s). The random distribution allowed a certain amount of duplication (to test whether a good candidate was consistently found with similar starting periods) to avoid grid-related issues. To evaluate a rough significance of pulsation candidates, we used the probability to find a given number of $\sigma$ from a search over 1000 realizations of white noise, using the same GTIs as the original observations. No significant (>3$\sigma$, using the above criterion) new pulsations were found.

## 3.3. The “Jerk” Search and a New Detection in 2016

Finally, we took advantage of the recently developed “Jerk” search technique in accelsearch (Andersen & Ransom 2018). This time, since this technique uses both the first and second frequency derivatives, we used longer chunks of data, around 80 ks (when the observation was long enough). Thanks to the “Jerk” search, we found one new candidate pulsation in ObsID 90201037002.

**New orbital solution**—We refined the candidate from the “Jerk” search with prepfold, and local values of the period and two period derivatives were measured with precision. They were consistent with those expected from an ~1.37 s pulsation Doppler-shifted by orbital motion with $P_{\text{orb}} \sim 2.52$ days and a $\sin i/c \sim 22$ lt-s, the orbital parameters known from B14. We refined the spin solution using the full observation (instead of an 80 ks segment) using HENphaseogram.<fnref n="24" /> These local spin derivatives give effects 2 orders of magnitude above the reasonable interval for the spin-up of the source, which is $|\dot{v}| \lesssim 2 \times 10^{-10}$ Hz s$^{-1}$. Despite the observation covering less than $1/3$ of an orbit, the measured spin derivatives were sufficient to put a constraint on the orbital phase (Figure 1). We can expect intrinsic (torque-driven) spin-up or spin-down whose magnitude is up to some $10^{-10}$ Hz s$^{-1}$. The spin-up and orbital parameters are degenerate, and we include these considerations when calculating the improved error bar on the orbital period $P_{\text{orb}}$.

Nonetheless, the 2.5 yr lever arm between this new measurement and the original 2014 data set allows for the error bar on $P_{\text{orb}}$

<footnote n="21">²¹ Rather than calculating the $Z_n^2$ statistics on the single events, this software prefolds the data and then calculates the statistics using the phases of the profile with a weight given by the number of counts in each profile bin (Huppenkothen et al. 2019).</footnote>

<footnote n="22">²² We chose 200 s because it is approximately the error on $T_0$ that would produce a shift of ~half a pulsar rotation over an orbit.</footnote>

<footnote n="23">²³ For example, a small dead-time component that repeats with a frequency of ~890 Hz in data sets taken in CP mode, usually visible in very bright sources.</footnote>

<footnote n="24">²⁴ The signal-to-noise ratio in single intervals was not adequate to make more sophisticated timing using, e.g., PINT or tempo2 at this step.</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** (a) The orbital phase can be measured from the spin solution if a measurement of at least two spin frequency derivatives ($\dot{v}$, $\ddot{v}$) is available. Here we compare the best spin solution calculated by PRESTO (cyan, and blue in the time range covered by the data) and the expected $v$ and $\dot{v}$ change due to the orbital Doppler effect (red). At this phase of the orbit, the constraint on the ascending node passage is very accurate. (b) New detections in ObsIDs 90201037002 and 80002092004 using the all-data orbital solution in Table 2. The horizontal line indicates the $3\sigma$ detection level for the given number of trial frequencies.</figcaption></figure>

to be significantly reduced (Table 2). With this rough orbital solution, we were able to fold the events in subintervals of the new observation and calculate the arrival times of the pulses with the `fftfit` method (Taylor 1992) as implemented in `HENpha-seogram`. Then, we used PINT to refine the orbital and spin parameters, fitting for the frequency, first derivative, and orbital period. Given the small observation length, the spin derivative is not constrained by the fit; therefore, we fixed it at zero. Finally, we found a well-constrained solution and an improved value of $P_{\text{orb}}$. This new value is close to the one calculated from 2014 data, even if it is outside the range allowed by the quoted error bar on $P_{\text{orb}}$ from the solution by B14. We will discuss the possible implications in Section 4.3.

Later, we tested the possibility that an integer number of additional orbital periods could fit into the time range between ObsIDs 80002092011 and 90201037002. An error on the period of ~0.0001 day, as implied by an additional or missing full orbit in this time range, would produce a very obvious distortion of the 2014 solution, and this possibility can be safely neglected (see Figure 2, bottom panel).

**Spin-down**—The most intriguing and direct consequence of this new detection is that, contrary to all other PULXs found up to now, the source showed an average spin-down (Figure 4) of $\dot{v} \sim -5.8 \times 10^{-11}$ Hz s$^{-1}$ between 2014 and 2016. The spin-down observed in this time range is likely driven by the negative torque of an accretion disk. The alternating behavior between spin-up and spin-down suggests that the source is close to spin equilibrium. An additional component could be dipole radiation from a strong magnetic field. We discuss the implications of this finding in Section 4.3.

### 3.4. Deep Pulsation Searches with New Orbital Solution and Upper Limits

Based on the improved orbital solution described in Section 3.3, we ran a deeper search using the $Z_1^2$ (Rayleigh test) and $Z_2^2$ statistics (Buccheri et al. 1983). Given the unexpected spin-down, this new search included a much larger interval of frequencies but with a better constraint on the orbital phase. This time, we also ran simulations to evaluate the pulsed fraction upper limit for the nondetections. To do so, for each data set, we

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Top panel: residuals of the pulse TOAs for the 2014 ObsIDs 80002092004–09 using the models in Table 3. Black dots are the residuals from model 1 if we ignore the orbital motion, and blue dots include the orbital motion. Green dots add the $\dot{v}$ component (model 2) instead. Bottom panel: residuals of the pulsar TOAs using a timing model similar to model 3, with two spin derivatives (F0, F1, F2) but changing the orbital period to accommodate one more (purple) or one less (orange) full orbit between ObsIDs 80002092011 and 90201037002. A very obvious modulation at the orbital period, with increasing amplitude, appears in the residuals.</figcaption></figure>

<table><caption>**Table 2**
Updated Orbital Parameters of M82 X-2 Compared to the B14 Fit and a New Fit of the 2014-only Data</caption><thead><tr><th>Parameter</th><th>B14</th><th>2014-only Fit</th><th>All Data</th></tr></thead><tbody><tr><td>T<sub>0</sub> (MJD)</td><td>56,694.7327(1)</td><td>56,682.0661(3)</td><td>56,682.0661(3)</td></tr><tr><td>a sin i (lt-s)</td><td>22.225(4)</td><td>22.215(5)</td><td>22.215(5)</td></tr><tr><td>P<sub>orb</sub> (day)</td><td>2.53260(5)</td><td>2.53287(5)</td><td>2.532948(4)</td></tr></tbody></table>

**Note.** The quoted uncertainties represent approximate symmetrical $1\sigma$ confidence intervals as calculated by PINT. Values and error bars on $a \sin i$ and $T_0$ are fixed to the 2014-only value. The 2014-only quoted error bar is based on the orbital phase measured in ObsID 80002092011, with a 40 day lever arm with respect to $T_0 = \text{MJD} 56,682.0661$ and including the local first derivative in the fit. The full 2014–2016 solution is based on the orbital phase measured in Figure 1 and refined with PINT, fixing the first derivative to zero (see text).
<table><caption>**Table 3**
Spin Parameters for the Observations in 2014 Using Different Models</caption><thead><tr><th>Parameter</th><th>Model 1</th><th>Model 2</th><th>Model 3</th><th>Model 4</th></tr></thead><tbody><tr><td>F0 (Hz)</td><td>0.728566160(7)</td><td>0.72852171(3)</td><td>0.72855082(8)</td><td>0.72854408(11)</td></tr><tr><td>F1 (10<sup>-11</sup> Hz s<sup>-1</sup>)</td><td>0</td><td>4.758(3)</td><td>-2.98(2)</td><td>-5.29(3)</td></tr><tr><td>F2 (10<sup>-17</sup> Hz s<sup>-12</sup>)</td><td>0</td><td>0</td><td>8.66(2)</td><td>10.93(3)</td></tr><tr><td>GL_F0 (Hz)</td><td>0</td><td>0</td><td>0</td><td>1.75(2) &times; 10<sup>-5</sup></td></tr><tr><td>GL_EP (MJD)</td><td>...</td><td>...</td><td>...</td><td>56,685.8 (fixed)</td></tr><tr><td>rms (&mu;s)</td><td>136</td><td>34</td><td>9.2</td><td>4.3</td></tr><tr><td>&chi;<sup>2</sup> (&mu;s<sup>2</sup>)</td><td>2698012</td><td>165709</td><td>12251</td><td>2704</td></tr></tbody></table>

**Note.** Orbital parameters were fixed to the values in Table 2. The reference epoch for all models is MJD 56,683. Note that, even if the F1 term is negative, suggesting spin-down, the effect of the second spin derivative F2 is dominant and gives spin-up.

maintained the GTIs and total number of photons in each interval and randomized the time of the events inside each GTI. There is no significant broadband noise at the frequency of the pulsar. Again, we used the statistical functions in *stingray* and *HENDRICS*. As discussed above, this software allows the statistics to be calculated from the folded profiles instead of the single events, as in the original Rayleigh test and Buccheri et al. (1983). To further speed up the calculation (in particular when calculating the larger number of realizations required by the upper limit calculations), we executed the total folding by folding *M* equal-length subintervals of the observations and aligning the folded subprofiles differently for slightly different values of the frequency (similar to the fast folding algorithm (Taylor & Weisberg 1982) and the technique used in *prepfold*). A bonus of this technique is that it is easy to shift the subintervals by integer values according to linear, quadratic, and higher-order laws and to measure rapid changes of trial spin frequencies. The main rules of thumb to avoid dispersion of the signal in multiple bins are that (1) the number of bins in the profile *N* is large enough that the pulse shape features are distinct in different bins, and (2) the number of profiles that are being shifted is at least twice the maximum shift of the subprofiles. We verified that calculating these statistical tests this way does not depart significantly from the expected statistics (more details in the *Appendix* and Figure 7). Note that this new folding algorithm, as of 2019 August 13, was made available in the development version of *HENDRICS*.

Using this method, we do not find significant detections in the new observations. As expected from the decrease of trial values and the better description of the data given by the orbital solution instead of a few spin derivatives, the detection in 90201037002 turns out to be highly significant (Figure 1). We only find a hint of pulsations in ObsID 30101045002, close to the predicted spin-down line and the detection limit, that we plot in Figure 4 with a gray point.

Surprisingly, despite not finding the pulsations in any more new observations, we did find a new significant detection in ObsID 80002092004, the second observation of the 2014 campaign, where pulsations were not detected by B14. The reason this happened is probably because the pulsation departs strongly from the smooth timing noise that was measured back then, shortly after the start of ObsID 80002092006, as we will discuss in Section 3.5.

We also repeated the search after filtering the events more aggressively by energy, only considering photons in the 8–30 keV interval. While in ObsID 8000209200X, M82 X-2 dominated the total flux, and all emission between 3 and

30 keV was significantly pulsed, the pulsations in ObsID 90201037002 are only detectable above 8 keV (which is understandable by looking at the spectra reported by Brightman et al. 2016 and following papers). This new search, besides improving the marginal detection in ObsID 30101045002, did not produce more detections. The full summary of detections and upper limits is reported in Figure 4.

### 3.5. A *Glitch*

Assuming the orbital solution given in Table 2, we fit the data with timing models including an increasing number of spin derivatives (Table 3). After adding the F2 ($\dot{v}$) term, besides additional long-term trends that might, in principle, be described by additional spin derivatives, an abrupt spin frequency change becomes apparent between ObsID 80002092004 and the start of 80002092006. Using the *glitch* plug-in of TEMPO2, we calculate two local spin derivatives (F0 and F1) from the pulse times of arrival (TOAs). Figure 3 shows the result: we find that the timing solutions need either a sudden frequency jump or a very strong frequency derivative at around the same time. The measured frequency change is too sharp to be described by a small number of additional spin derivatives. It amounts to around $10^{-5}$ Hz, corresponding to $\Delta\nu/\nu \approx 1.3(4)10^{-5}$ (Table 3). The minimum $\dot{v}$ necessary to produce such a large frequency increase on timescales of 1 day is $\sim 10^{-10}$ Hz s$^{-1}$. Based on simple arguments from accretion theory, such a strong frequency derivative would need a significant change of mass accretion rate and, presumably, luminosity that we do not observe. We therefore discuss from now on the simplest explanation, i.e., that the neutron star has undergone a “glitch,” a sudden change of spin velocity.

The measured frequency change of $\Delta\nu/\nu \gtrsim 10^{-5}$ is a very large value, but not unprecedented, for a pulsar glitch. Typical glitch magnitudes in radio pulsars are more than an order of magnitude (often many orders of magnitude) smaller, but those observed in magnetars and accreting pulsars can reach those values; see Manchester (2018) for a review and Section 4.4 for a discussion.

To exclude instrumental effects, we checked the details of the clock correction file during this time range. We compared the results of the standard clock correction pipeline with a new clock correction based on an improved thermal model of the spacecraft oscillator (M. Bachetti & C. B. Markwardt 2020, in preparation), finding a very good agreement (well inside the error bars). Bad clock offset measurements from the ground stations could, in principle, produce changes in the measured
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Top panel: TOA phase residuals after correcting the timing solution with model 3 in Table 3. Middle panel: local frequency variation calculated with the *glitch* tool of tempo2 using nonoverlapping subsets of TOA residuals. The rapid change of frequency around MJD 56,685.7 (vertical line), ~10 µHz, is clearly off with respect to the much smoother (but implying already quite extreme spin-up) frequency variation at later times. Bottom panel: *v* calculated from the values of *v* in the middle panel. The minimum average *v* necessary to produce the observed change of frequency between MJD 56,685 and 56,686.5 is ~8 × 10⁻¹¹ Hz s⁻¹, very much above the local values of *v*.</figcaption></figure>

pulsed frequency, but we verified that the scatter of the offset measurements from the Malindi station is ≤100 µs, and we excluded all measurements from other, less reliable ground stations with no significant changes in the results.

## 4. Discussion

The strong spin-down observed between 2014 and 2016, the dramatic spin-up and pulsed fraction variation during the 2014 observation, and the nondetection of pulsations in subsequent observations, even with comparable fluxes of M82 X-2 measured by *Chandra*, can shed light on M82 X-2 and the ULX phenomenon on the whole. The glitch revealed in the 2014 observation also adds to the list of peculiar phenomena observed in ULXs, following the detection of antiglitches in NGC 300 X-1 (Ray et al. 2018).

In this section, we discuss how these new data can shed light on the nature of M82 X-2 and the ULX phenomenon.

### 4.1. Caution: Luminosity and Accretion Rate in M82 X-2

Most of the interpretation in the following subsections, in particular for what concerns the transient pulsations and the spin-up, could be done by simply assuming that the luminosity is a simple function of mass accretion rate.

However, this needs a few words of caution. An example of why this simple assumption might be flawed can be found in Vasilopoulos et al. (2019): they found that, despite a large drop of luminosity of the source NGC 300 ULX-1, the pulsar keeps spinning up with the same rate as when the luminosity is high.

This is probably due to the fact that the luminosity drop is given by occulting material, like a super-Eddington wind, and not an actual change of accretion rate.

If we take for granted that the drop of flux of M82 X-2 corresponds to a change of accretion rate, the finding by Brightman et al. (2019) that this flux change is periodic would imply that the mass accretion rate changes periodically, switching accretion on and off every ~60 days. For the physical properties of this system, however, this is problematic. While it would be easy to produce a periodically changing mass transfer rate if the period corresponded to the orbital period, this is far from obvious for a superorbital periodicity. Assuming Roche lobe overflow, possible mechanisms to periodically increase the pressure at the L1 Lagrangian point might involve additional orbiting objects, like the so-called Kozai mechanism (Kozai 1962; Lidov 1962) or tidally excited oscillations (e.g., Fuller 2017). Note that a very small periodic change of orbital eccentricity, well below the quoted upper limit from B14, would be sufficient to produce substantial changes of accretion rate (e.g., Hut & Paczynski 1984; Maccarone et al. 2010). In the same way, a stellar pulsation could produce a periodic increase of the mass accretion rate. However, we can exclude an evolved star due to the small semimajor axis of the orbit (Fragos et al. 2015; Heida et al. 2019) and the likely range of companion masses derived from timing (5–20 $M_\odot$; B14). For this range of masses, we expect the companion to be a main-sequence O-type star, and pulsations reported for these systems (e.g., from β Cepheids or slowly pulsating Be stars) are on shorter timescales. (see Samus’ et al. 2017). Stars that can, in principle, have stable pulsations in the 60 day range are usually either evolved lower-mass stars or luminous blue variables (LBVs) that are thought to be much more massive than the likely 5–20 $M_\odot$ range for M82 X-2 (see, e.g., Conroy et al. 2018; Jiang et al. 2018). Note also that Heida et al. (2019) ruled out LBVs in the error circle of M82 X-2.

The point of these arguments is that historic luminosity changes in M82 X-2 might be dominated by some occulting material, such as a precessing disk, and not by actual mass accretion rate changes. The variability of M82 X-1 complicates the investigation further and undermines a precise interpretation of the pulsar behavior based on luminosity.

Future studies will address this problem. For now, we keep both doors open and discuss the consequences of both hypotheses.

### 4.2. Pulsed Fraction Evolution and Nondetections

The pulsed fraction of M82 X-2 changed dramatically in 2014, and in the new observations, it often showed no detectable pulsations. In Table 1, we report the estimated pulsed fraction and/or upper limits on pulsations, plus the luminosity from the PULX M82 X-2 and the flux ratio between M82 X-2 and M82 X-1 when *Chandra* observations are available. In the latter cases, we also correct the upper limit of the pulsed fraction. The detectability of pulsations is largely driven by the intensity of the source *I*, the background *B*, and the pulsed fraction or, equivalently, the rms *r*, and only weakly by the observing time *T*. Lewin et al. (1988) showed that the significance of pulsations goes roughly as $I^2/(I+B)r^2T^{1/2}$. Looking in detail at the table, we see that the 2014 observations were executed in a very favorable situation, with M82 X-2
<figure><img src="image_4.png" /><figcaption>**Figure 4.** History of detections of pulsations from M82 X-2, showing the mean spin-down of the neutron star. The new detections are highlighted with larger symbols (including the marginal one in gray). The inset shows the timing solution from model 3 in Table 3, applied to the data from ObsIDs 80002092004–09. Upper limits are indicated by arrows. Note that for the marginal detection in 30101045002, we plot both the marker and the upper limit.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Each point represents 80 ks of data during the 2014 campaign. The pulsed fraction increases with total flux. This implies that the total flux increase is driven by M82 X-2, not M82 X-1, and allows us to assume intensity as a proxy for accretion rate for the 2014 data. Also, the increase of the total pulsed fraction (pulsed flux over the total flux of M82 X-1 and M82 X-2) is too steep to be justified by just the flux evolution of M82 X-2 over M82 X-1 (the dashed lines show the expected increase in this case). There seems to be an approximately linear increase of pulsed fraction with M82 X-2's intensity (blue line). Hardness increases slightly with intensity, as expected from an increase of M82 X-2, which is generally harder (Brightman et al. 2016).</figcaption></figure>

brighter than M82 X-1 during ObsID 80002092007, a high (even if not the highest value) flux of M82 X-2 and a very long observing time. The pulsed fraction change in 2014 is positively correlated with a flux increase of M82 X-2 (see below), and we could have expected even higher values in these later observations.

Figure 5 shows that the increase of the pulsed fraction is associated with an increase of flux and a stable hardness ratio. This tells us that the pulsed fraction is increasing because M82 X-2 is brightening and not, e.g., because of a weakening of M82 X-1. Therefore, it is safe to associate the change of

luminosity with a change of accretion rate and interpret the increasing torque implied by the strong second derivative $\dot{v}$ in these terms as we do in Section 4.3.

On the contrary, a few observations in the following years were clearly in unfavorable conditions for pulsation detection. In ObsIDs 80202020002, 30202022002, and 90202038004, the flux of M82 X-2 was extremely low. In 2015–2016, M82 X-1 underwent some flaring events, and pulsations would have been impossible to detect in ObsIDs 30202022008 and 90202038002 despite the source luminosity being $\gtrsim 10^{39}$ erg s$^{-1}$. The low pulsed fraction measured in ObsID 90201037002 probably implies that M82 X-1 was strong with respect to M82 X-2 but not strong enough to completely drown the pulsations. On the other hand, the detection of pulsations tells us that M82 X-2 was stronger than in the two observations performed approximately 2 weeks before and after, 30202022010 and 90202038002, despite the overall NuSTAR count rate being similar.$^{25}$

In ObsIDs 90101005002, 80202020004, and 80202020008, the flux of M82 X-2 was higher than in 2014. The flux of M82 X-1 was also higher in these observations, and we would not expect the same significance as 2014, but still, we did not find credible candidate pulsations, even marginal, in these ObsIDs, even if the pulsed fraction upper limit was compatible with the values measured in 2014. Also, 80002092002 and the start of 80002092004 show no detectable pulsations. This might be due to the onset of accretion, with the higher flux in 80002092002 related mostly to M82 X-1. The ULX M82 X-2 is not the only one showing intermittent pulsations. Recently, Sathyaprakash et al. (2019) found intermittent pulsations in NGC 1313 X-2. The ULXs are peculiar objects, probably Roche lobe-filling like low-mass X-ray binaries (that have much lower companion masses) but with companion masses

<footnote>$^{25}$ Note that the change of flux from M82 X-2 on these timescales is compatible with the superorbital period of 60 days shown by Brightman et al. (2019).</footnote>
more similar to high-mass X-ray binaries (that mostly accrete via wind).<fnref n="26" />

Therefore, it is useful to look for transient pulsations in both kinds of systems and see what processes might be at work. One notable example is Aql X-1, an accreting millisecond pulsar in a low-mass X-ray binary, probably with a relatively low magnetic field (10⁹ G) and sub-Eddington. Casella et al. (2008) reported pulsations lasting only ~150 s over several megaseconds of RXTE observations. Another accreting millisecond pulsar, HETE J1900.1-2455, showed transient pulsations during an outburst, with changes of pulsed fraction associated with thermonuclear bursts (Galloway et al. 2006). The HMXB A0538-66 showed pulsations in 1982 (Skinner et al. 1982) but none since (see Kretschmar et al. 1997, for a review). Recently, Brumback et al. (2018) and Pike et al. (2019) reported on transient pulsations from the X-ray binaries LMC X-4 and SMC X-1. In the first case, a significant change of pulse properties (rms, spin derivative) was associated with the precursor of a large burst-like flux increase, which is compatible with an increase of the accretion rate at the surface (e.g., the disk overcoming the centrifugal barrier of the magnetized neutron star). In the second, there was, instead, evidence of a change of absorption between the nonpulsed and pulsed intervals, suggesting the presence of occulting material.

A precessing structure that occults the central X-ray source could be the explanation for the periodic modulation of the flux from M82 X-2 (Middleton et al. 2015b; Brightman et al. 2019). A strong disk wind, as observed in several ULXs (e.g., Fabrika et al. 2015; Pinto et al. 2017; Kosec et al. 2018), would follow a precessing disk (e.g., Pasham & Strohmayer 2013; Middleton et al. 2019b) and block the view of the central flux from the compact object. For example, Dauser et al. (2017) modeled the long superorbital periods observed in ULXs as the precession from a disk launching a conical wind. Precession is observed in an archetypal supercritical source, SS 433, and various ULXs (Begelman et al. 2006). It is expected to produce a periodic hardening of ULX spectra (Middleton et al. 2015a).

Using the best Chandra data available, Brightman et al. (2020), found no significant spectral evolution of M82 X-2 and, in particular, no evidence of large variations of $N_H$. The spectrum is often consistent with the pulsed spectrum reported by Brightman et al. (2016). However, the $N_H$ is quite high, and the spectrum of M82 X-2 is hard, with the Galactic X-ray emission dominating the lowest energies below the NuSTAR band. It is possible that changes of $N_H$ go unnoticed or confused with other degenerate parameters. Also, a highly ionized wind would reduce the flux by scattering rather than the photoelectric effect, reducing the measured $N_H$.

An additional explanation might involve the precession of the pulsar, with the pulse beam moving away from the observer. In principle, this might be detected through changes of the pulsar spectrum. Given the relatively low count rate of the source and the large source confusion (with M82 X-1 at 5"), this is not testable on NuSTAR data and difficult even from existing Chandra spectral data that are sparse and often affected by pileup, which affects the higher energies where the pulsed component is more significant. In principle, the large spin-down discussed in the next section might indicate a change of the accretion rate that could indeed alternate between

<figure><img src="image_6.png" /><figcaption>Figure 6. Using the same data as Figure 5, we compare the measured values of intensity and spin-up with the predictions of torque theories in the case of a slow rotator, where $R_{in} \ll R_{co}$. Far from spin equilibrium, we would expect the frequency derivative to be proportional to $\dot{M}^{6/7} \propto L_{37}^{6/7}$ (Frank et al. 2002). For slim disks, $L \propto \dot{M}^2$ (King 2008; King et al. 2017). The data are incompatible with the slow rotator hypothesis, suggesting that $R_{in} \sim R_{co}$. Indeed, using Equation (4) with $\omega_{fast} \sim 1$, the curve can be made to agree with the data.</figcaption></figure>

accretion and no accretion, with the total emission dominated by a disk close to spin equilibrium. With the available data, it is not possible to get a definitive answer.

### 4.3. Spin Behavior: Not a Slow Rotator

The absence of pulsations for a high fraction of time and the secular spin-down observed between 2014 and 2016, despite luminosity trends consistent with past observations, suggest that the source might be close to spin equilibrium, with large variations of spin-up or spin-down produced by relatively small changes of the magnetospheric radius $R_M$ around the corotation radius $R_{co}$.

Using the formulae by Ghosh & Lamb (1979b) and assuming constant $\dot{M}$, one can show that the spin-up of a pulsar should follow the relation

$$-\frac{\dot{p}}{p} \propto p n(\omega_{\text{fast}}), \qquad (3)$$

where $p$ is the pulse period, $\dot{p}$ is its first derivative, $n$ is the torque from the disk on the magnetic field lines, and $\omega_{fast} = (R_M/R_{co})^{3/2}$ is the fastness parameter. Wang (1995) derived the following expression for the torque:

$$n(\omega_{\text{fast}}) \approx \frac{7/6 - (4/3)\omega_{\text{fast}} + (1/9)\omega_{\text{fast}}^2}{1 - \omega_{\text{fast}}}. \qquad (4)$$

If the pulsar is a slow rotator, with $\omega_{fast} \ll 1$ and thus very far from spin equilibrium, this relation reduces to a constant, $n(\omega_{fast}) \approx 7/6$, and Equation (3) gives

$$-\frac{\dot{p}}{p} \propto p. \qquad (5)$$

This relation holds, for example, in NGC 300 ULX-1 (Vasilopoulos et al. 2018, 2019).

However, the spin-up shown by M82 X-2 during the 2014 campaign does not follow this relationship at all (Figure 6). If we consider a change of $\dot{M}$ proportional to the change of luminosity of the source, we expect a relation of the kind (see

<footnote n="26">²⁶ For completeness, we mention that Mellah et al. (2019) proposed an intermediate Roche wind accretion model, initially proposed by Mohamed & Podsiadlowski (2007), for ULXs.</footnote>
Ghosh & Lamb 1979b)

$$-\frac{\dot{p}}{p^2} = \dot{v} \propto L_{37}^{6/7} n(\omega_{\text{fast}}), \qquad (6)$$

while in the “slim disk” case with beaming (King 2008, 2009), $\dot{M} \propto L^{3/7}$ and

$$-\frac{\dot{p}}{p^2} = \dot{v} \propto L_{37}^{3/7} n(\omega_{\text{fast}}). \qquad (7)$$

Again, Figure 6 clearly shows for both models that the spin behavior of M82 X-2 in 2014 is inconsistent with $\omega_{\text{fast}} \ll 1$, while it can, in principle, be reproduced if we assume $\omega_{\text{fast}} \sim 1$. Therefore, our results support the basic assumption of most concurring models on M82 X-2 (e.g., B14; Dall’Osso et al. 2015; Tsygankov et al. 2016; King et al. 2017). Here we calculate spin derivatives from a robust interpolation<fnref n="27" /> of order 2 of single spin frequency measurements during the campaign.

The source might be undergoing some sort of spin reversal over time. This would be similar to what was reported for LMC X-4 by Molkov et al. (2017), where an alternating spin-up and spin-down behavior was observed over many decades. In that case, one of the possible explanations was a recycling magnetosphere model (Perna et al. 2006), where accretion is only possible at certain spin phases.

Independent of the processes investigated above, the pulsar might be spinning down because of its own rotation-powered radiation. Given this source’s large spin-up in other observations, this mechanism is unlikely to dominate over accretion torques; however, for the sake of completeness, it is interesting to make an estimate of the magnetic field in this case. According to the classical magnetic dipole radiation formula,<fnref n="28" /> if the spin-down of the pulsar is given by dipole radiation, we can estimate the magnetic field as

$$B > \left( \frac{3c^3I}{8\pi^2R^6} \dot{p}\dot{v} \right)^{1/2} \approx 3.2 \times 10^{19} \left( \frac{\dot{p}\dot{v}}{s} \right)^{1/2} \text{G}, \qquad (8)$$

where I is the moment of inertia of the neutron star and R is its radius. Substituting the spin period and spin-down period derivative in the formula above and using standard estimates for the radius and moment of inertia of the neutron star, we get a lower limit on the magnetic field of $\sim 3 \times 10^{14}$ G if the dipole axis is perpendicular to the rotation axis. Again, this assumes that there is no torque-driven spin-down and that the emission is from a “classical” rotating misaligned magnet. The actual radio pulsar emission process is still unknown (see Pétri 2016, for a review) and is believed to involve collective plasma motions that might be difficult to produce in an accreting system.

Future robust measurements of the total mass exchange will give a better estimate of the position of the spherization radius and help in understanding the viability of alternative explanations to the accretion-driven spin-down.

## 4.4. Glitch: Strong Spin-down before 2014?

There are very few glitches reported from accreting neutron stars. In the early 2000s, KS 1947+300 showed a large glitch of $\Delta\nu/\nu \approx 3.5 \times 10^{-5}$ in RXTE observations (Galloway et al. 2004), comparable with $\Delta\nu/\nu \approx 10^{-5}$ reported in Section 3.5. Similarly, Serim et al. (2017) reported on a strong glitch from the slow accreting pulsar SXP 1062, with $\Delta\nu/\nu \approx 10^{-3}$. Besides SXP 1062 and KS 1947+300, glitches with $\Delta\nu/\nu \gtrsim 10^{-5}$ have only been reported in magnetars (e.g., Dib et al. 2009; Dib & Kaspi 2014; Archibald et al. 2016). Some sources, notably the ultraluminous pulsar NGC 300 ULX-1 (Ray et al. 2018), show antiglitches instead—sudden slow-downs of the pulsar. A nonmagnetar glitch is believed to be due to the independent rotation of the neutron star core with respect to the rest of the neutron star. A neutron star contains superfluid material in the core and inner crust, and the different regions of the star interior are allowed to rotate with different velocities. If the crust slows down over time, e.g., due to dipole radiation as in radio pulsars, the spin rates of the crust and core differ more and more over time. It is believed that when the differential rotation reaches some critical value, these regions connect, angular momentum is transferred between the crust and the core, and the star suddenly starts spinning as a whole. If the crust was slowing down over time, the connection with the core will produce a sudden acceleration. The opposite is believed to have happened with the antiglitch observed in NGC 300 ULX-1 (Ray et al. 2018). In that case, the star has been observed to spin-up over time due to accretion, and the antiglitches likely represent the effect of reconnection between the faster crust and slower core. It is possible that the glitch we observe at epoch 56,685.8 is a standard glitch. This would require that the star was spinning down prior to the 2014 observations. Given the spin-down reported between 2014 and 2016, it is indeed possible that a similar spin-down was acting on the neutron star before that time, leading to the glitch. Another option would be that the process involved here is due to a high magnetic field, like in magnetars. This would agree with the secular spin-down being due to dipole radiation. We do not have enough data to test this hypothesis.

## 4.5. An Orbital Period Derivative?

As reported in Section 3.3 and Table 2, if we select only the data from the 2014 observations and fit the orbital parameters, we obtain a slightly different orbital period than measured, including the 2016 data set. This might, in principle, mean that the orbit has shrunk and the orbital period has decreased over time, due, for example, to the strong mass exchange that the super-Eddington luminosity of the source seems to imply. However, the 2014 observations contain a strong timing noise that might have influenced the fit of the orbit reported by B14 and adjusted in this paper. We find slightly different solutions depending on the model for the red noise, and Table 2 sums up this variability through larger error bars than the ones calculated from the single fits.

However, if the reported tension between the 2014-only orbital solution and the one with the full data set is true, a single new measurement of the orbital phase in future observations (at this point, more than 3 yr since the last one) will be able to measure with high significance a negative orbital period derivative. The values in Table 2 would imply $\dot{P}_{\text{orb}} \sim -10^{-7}$ days day$^{-1}$ implying an orbital evolution timescale of 30,000 yr and a mass exchange almost 3 orders of magnitude above the Eddington limit (see below).

<footnote n="27">²⁷ We use the method by Savitzky & Golay (1964), as implemented in scipy.signal.</footnote>

<footnote n="28">²⁸ See any textbook on pulsars, e.g., https://www.cv.nrao.edu/course/astr534/Pulsars.html.</footnote>
Following Rappaport et al. (1982), the expected orbital period derivative for mass transfer from a more massive companion star, neglecting a contribution from gravitational waves, can be estimated through

$$ \frac{\dot{a}}{a} = \frac{2\dot{P}_{\text{orb}}}{3P_{\text{orb}}} = -\frac{\dot{M}_{c}}{M_{c}}\left[2\left(1-\frac{\beta}{q}\right)-\frac{1-\beta}{1+q}-\frac{2\alpha(1-\beta)(1+q)}{q}\right], \quad (9) $$

where $M_c$ and $\dot{M}_c$ are the mass and mass-loss rate from the companion (negative); $M_p$ and $\dot{M}_p$ are the same quantities for the pulsar, $q = M_p/M_c$, $\beta = -\dot{M}_p/\dot{M}_c$; and $\alpha$ is the specific angular momentum of the lost mass in units of $2\pi a^2/P_{\text{orb}}$. In the conservative scenario, $\beta = 1$, and Equation (9) becomes

$$ \frac{\dot{a}}{a} = \frac{2\dot{P}_{\text{orb}}}{3P_{\text{orb}}} = -\frac{\dot{M}_c}{M_c}\left[2\left(1-\frac{1}{q}\right)\right]. \quad (10) $$

In the scenario where all mass is lost in an outflow before accreting onto the pulsar, but the outflow is launched from close to the pulsar, the specific angular moment is $j_{\text{outflow}} \approx 2\pi a_p^2/P_{\text{orb}}$, where $a_p$ is the semimajor axis of the pulsar orbit. Therefore, $\beta = 0$ and

$$ \alpha \approx \left(\frac{a_p}{a}\right)^2 = \left(\frac{M_c}{M_p + M_c}\right)^2 = \frac{1}{(1+q)^2}, $$

and Equation (9) reduces to

$$ \frac{\dot{a}}{a} = \frac{2\dot{P}_{\text{orb}}}{3P_{\text{orb}}} = -\frac{\dot{M}_c}{M_c}\left[2 - \frac{2+q}{q(1+q)}\right], \quad (11) $$

which is a relatively small correction to Equation (10) for small values of $q$.

If the mass transfer in the system is highly super-Eddington and no significant mass losses happen close to the donor star, both scenarios predict that the orbit shrinks, with orbital period derivatives of the order

$$ \dot{P}_{\text{orb}} \sim -0.015 \left( \frac{M_p}{1.4 M_{\odot}} \right)^{-1} \left( \frac{-\dot{M}_c}{M_{\odot}/\text{yr}^{-1}} \right) \text{s s}^{-1} \quad (12) $$

$$ \sim -3.5 \times 10^{-8} \left( \frac{M_p}{1.4 M_{\odot}} \right)^{-1} \left( \frac{-\dot{M}_c}{100 \dot{M}_{\text{Edd}}} \right) \text{s s}^{-1}. \quad (13) $$

This $\dot{P}_{\text{orb}}$ corresponds to a change of orbital period of $\sim 1 \text{s yr}^{-1}$ and produces a shift of the orbital phase by a few hundred seconds over a few years. Again, since $\dot{M}_c$ is negative, this $\dot{P}_{\text{orb}}$ is also negative. This is easy to observe with standard techniques of pulsar timing. Future observations will test this hypothesis by tracking the orbital phase over time.

Finally, in the remaining extreme scenario, where most mass is lost from the system in the form of winds from the donor star, we have

$$ \frac{\dot{a}}{a} = \frac{2\dot{P}_{\text{orb}}}{3P_{\text{orb}}} = -\frac{\dot{M}_c}{M_c}\frac{1}{(1+q)}. \quad (14) $$

It is easy to show that this scenario would cause the orbit to expand instead of shrink (positive $\dot{P}_{\text{orb}}$).

It is likely that in a real-life scenario, both phenomena should happen, stabilizing the mass transfer over long timescales. Fragos et al. (2015) discussed these arguments in detail. The detection of an orbital period derivative would set strong constraints on the viable models for the evolution of the binary system in M82 X-2.

## 5. Conclusions

This work characterizes the timing behavior of M82 X-2 over 2 yr, showing a number of phenomena that can be used to understand the nature of this remarkable source. None of these new findings provide definitive information on the nature of the source, but they represent important clues. Thanks to this work, we now know the following.

1. The pulsar alternates phases with a large spin-up with phases with a large spin-down. Moreover, the spin-up is inconsistent with small values of the fastness parameter. This is a strong indication that the pulsar is close to spin equilibrium.

2. The pulsations are transient, changing their significance over time, with a tentative positive correlation with the flux of the source. This points to an intrinsic mechanism for this pulsed fraction variability rather than, e.g., a change of the background flux from the nearby M82 X-1.

3. The neutron star has very strong glitches, probably due to the rapid spin evolution. The observed positive glitch suggests that a strong spin-down occurred prior to the first detection in 2014.

In addition, we now have a very precise orbital solution that can be used to look for orbital phase derivatives in future observations. This will be crucial to test the total mass exchange in the system.

M.B. thanks the Fulbright Visiting Scholar Program for supporting a 9 month visit at Caltech. M.M. appreciates support from an STFC Ernest Rutherford Fellowship. Part of the pulsar search and optimization software was developed in the framework of the Citizen Computing Pulsar Search (CICLOPS), a project supported by POR FESR Sardegna 20142020 Asse 1 Azione 1.1.3 (code RICERCA_1C-181), call for proposal “Aiuti per Progetti di Ricerca e Sviluppo 2017,” managed by Sardegna Ricerche. The authors wish to thank Jim Fuller for insightful conversations on stellar pulsations, Georgios Vasilopoulos for comments on the manuscript, Sergei Popov for pointing out alternative models for magnetar evolution, and Wynn Ho for pointing out an incorrect citation about glitch models in the first version posted on the arXiv. Thanks also to Deepto Chakrabarty, Fred Lamb, Juri Poutanen, Alexander Mushtukov, and Andrew King for multiple insights on accretion theories and competing models for PULXs.

*Facilities: NuSTAR* (Harrison et al. 2013), *Chandra* (Weisskopf et al. 2002), ATNF pulsar catalog (Manchester et al. 2005).

*Software: HEAsoft (HEASARC 2014), FTOOLS (Blackburn et al. 1999; Blackburn 1995), Stingray (Huppenkothen et al. 2016, 2019), HENDRICS (Bachetti 2018), astropy (Astropy Collaboration et al. 2018), PINT (Luo et al. 2019), PRESTO (Ransom 2011).
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Statistical distribution of $Z_1^2$ values in our blind search of pulsations, using the orbital solution in Table 2 and a search over $T_0$ values and $f, \dot{f}$. The dotted, dashed-dotted, and dashed vertical lines show the $3\sigma$, $4\sigma$, and $5\sigma$ levels from the theoretical distribution for the case of a single search over $f - \dot{f}$ and accounting for the given number of trial values of $T_0$. Panels (a) and (b) show ObsID 80002092004 in 2014. The detection significance is barely above $3\sigma$ if considered as a blind search, but the spin values and inferred $T_0$ are very close to the values in ObsID 80002092006, which increases the reliability of this detection. Panels (c) and (d) show the same plot for the new detection in ObsID 90201037002, showing a detection significance well above $5\sigma$.</figcaption></figure>

## Appendix

### On the Detection Limits Used in the Pulsation Search

Figure 7 shows the statistical distribution of $Z_n^2$ values during a search for pulsations with the methods described in Section 3.4. We executed the analysis on real and (hundreds of realizations of) simulated data with the same GTIs and count rate as the original observation. In a search containing only simulated white noise, the statistical distribution very closely follows the expected $\chi^2$ distribution (Buccheri et al. 1983). This means that the folding-and-shift procedure, with a careful choice of parameters, does not significantly alter the null hypothesis probability that can be used to search for pulsations.

We conclude that the slight deviations from the white-noise statistical distribution found in real data are driven by the source variability and not the analysis technique used.

To evaluate the p-value, in order not to synthetically increase the number of trials, we consider only the trials over $T_0$ needed to get the best detection: if the best $T_0$ is very close to the known orbit, it means that the search over the full $T_0$ space is noise without any added value.

### A.1. Parameter Files

In this section, we write down a few useful parameter files to show the logical steps to obtain the final orbital solution.

The first step consists of fitting the two ObsIDs with the most significant pulsations (80002092007 and 80002092009) and covering multiple orbits, using a spin solution with four spin derivatives, in order to constrain the orbital period and, in particular, the ascending node passage. The result is the following set of timing parameters (in TEMPO2/PINT format):

<table><tr><td>PSR</td><td>M82_X2_0BSID_007_009</td></tr><tr><td>RAJ</td><td>9:55:42.14400000</td></tr><tr><td>DECJ</td><td>69:40:26.00400000</td></tr><tr><td>F0</td><td>0.72856670954598399391<br>8.741567555166416003e-08</td></tr><tr><td>F1</td><td>6.4250081648893346535e-111<br>4.9287115294950599775e-13</td></tr><tr><td>F2</td><td>2.2525682844363675388e-161<br>7.8038412567833659266e-18</td></tr><tr><td>F3</td><td>-5.939687798010528739e-221<br>5.0019855875541822192e-23</td></tr></table>
<table><thead><tr><th></th><th>(Continued)</th></tr></thead><tbody><tr><td>PEPOCH</td><td>56695.507895798600000</td></tr><tr><td>PLANET_SHAPIRO</td><td>N</td></tr><tr><td>BINARY BT</td><td></td></tr><tr><td>PB</td><td>2.5328840644740456124 1<br>0.00016109793561046761965</td></tr><tr><td>A1</td><td>22.204213222811447 1<br>0.012032723626796296</td></tr><tr><td>ECC</td><td>0.0</td></tr><tr><td>T0</td><td>56682.068219262420266 1<br>0.001018781929672351325</td></tr><tr><td>OM</td><td>0.0.</td></tr></tbody></table>

Once we do that, we can use ObsID 80002092011, fixing $T_0$ and using the longer lever arm to better estimate the orbital period $P_{orb}$:

<table><tbody><tr><td>PSR</td><td>M82_X2 OBSID_011</td></tr><tr><td>RAJ</td><td>9:55:42.14400000</td></tr><tr><td>DECJ</td><td>69:40:26.00400000</td></tr><tr><td>F0</td><td>0.72875923701101958546 1<br>3.364666714127084281e-07</td></tr><tr><td>F1</td><td>1.189117392632877145e-10 1<br>6.9989856290673740436e-12</td></tr><tr><td>PEPOCH</td><td>56720.878600229010000</td></tr><tr><td>PLANET_SHAPIRO</td><td>N</td></tr><tr><td>BINARY BT</td><td></td></tr><tr><td>PB</td><td>2.5328387434901325451 1<br>4.4959856101332736067e-05</td></tr><tr><td>A1</td><td>22.204213222811447 0<br>0.012032723626796296</td></tr><tr><td>ECC</td><td>0.0</td></tr><tr><td>T0</td><td>56682.068219262420266 0<br>0.001018781929672351325</td></tr><tr><td>OM</td><td>0.0.</td></tr></tbody></table>

This is the solution reported as “2014-only” in Table 2.

Finally, we can use the 2010 ObsID 90201037002 to further constrain the orbital period, obtaining

<table><tbody><tr><td>PSR</td><td>M82_X2 OBSID_90201037</td></tr><tr><td>RAJ</td><td>9:55:50.20800000</td></tr><tr><td>DECJ</td><td>69:40:46.99200000</td></tr><tr><td>F0</td><td>0.72390599836134799007 1<br>7.123714905580943099e-07</td></tr><tr><td>PEPOCH</td><td>57641.998517226645995</td></tr><tr><td>PLANET_SHAPIRO</td><td>N</td></tr><tr><td>BINARY BT</td><td></td></tr><tr><td>PB</td><td>2.5329477619141624702 1<br>4.0457956366063808576e-06</td></tr><tr><td>A1</td><td>22.204213222811447 0<br>0.012032723626796296</td></tr><tr><td>ECC</td><td>0.0</td></tr><tr><td>T0</td><td>56682.068219262420266 0<br>0.001018781929672351325</td></tr><tr><td>OM</td><td>0.0.</td></tr></tbody></table>

## A.2. Would a High Magnetic Field Imply an Electron-Capture Origin?

The findings of Sections 3.3 and 3.5 can be interpreted as signatures of a very high magnetic field. The idea of a magnetarscale magnetic field in a binary, however, is puzzling. A typical supernova explosion will result in the loss of a substantial fraction of the mass from the supernova progenitor, which will usually be

the heavier object in a binary. When more than half of the total mass in the system is lost, the binary becomes entirely unbound. With more modest mass loss, the binary becomes eccentric and acquires a longer orbital period than it had before the merger (e.g., Boersma 1961). In such systems, accretion can take place only for a small range of orbital phases very close to periastron, which would suppress accretion. Given that the typical decay timescale for a magnetar’s magnetic field is of order 10,000 yr,<fnref n="29" /> but the circularization timescale for donor stars with radiative envelopes is much longer than this, the combination of a nearly circular orbit with a magnetar-like magnetic field requires some explanation.

In principle, a very finely tuned asymmetry of the supernova explosion can “fix” some of the problems caused by mass loss, but these systems will have long orbital periods (Kalogera 1998). The alternative solution is to require a minimal mass loss. This can be done by electron-capture supernova processes, in which an iron core never forms, and can result either from the collapse of intermediate-mass stars (e.g., Nomoto 1987; Poelarends et al. 2008) or from the accretion-induced collapse of a white dwarf in a binary (Ergma 1993; Fryer et al. 1999). Some observations of accreting X-ray pulsars suggest that there are two sets of these objects, with reasonable—albeit not airtight—arguments that some form in core-collapse supernovae and others form in electron-capture supernovae (Pfahl et al. 2002; Knigge et al. 2011). These systems are expected to have not just lower symmetric kicks due to mass loss (Blaauw 1961) but also lower asymmetric kicks (Fryer et al. 1999).

A few other aspects of this scenario are additionally attractive. First, evolutionary calculations show a preference for donors of about $5–10 M_\odot$ (Fragos & McClintock 2015). The masses of stars in binaries tend to show correlations, with a flat mass ratio distribution, rather than the mass ratio distribution predicted by randomly drawing from a standard initial mass function law (Sana et al. 2012). The direct electron-capture supernovae tend to occur for progenitor masses of about $8–10 M_\odot$ (Nomoto 1987), while the white dwarfs with small enough carbon abundances that they can undergo accretion-induced collapses have progenitor masses just below the range required for supernovae to occur (Canal & Schatzman 1976). Furthermore, there are widespread suggestions that accretion-induced and/or merger-induced collapse events (i.e., those in which the merger of two white dwarfs drives a collapse to a neutron star) can lead to the production of magnetars (e.g., Usov 1992; King et al. 2001). Additionally, the lack of a supernova remnant associated with such a young system also favors the idea that the neutron star formed in a manner with much less mass ejection than typical core-collapse supernovae. If indeed the magnetic field is high, as suggested by the spin-down and glitch (and as previously suggested for this source in the literature, e.g., Dall’Osso et al. 2015; Ekşi et al. 2015; Mushtukov et al. 2015; Tsygankov et al. 2016), there is, thus, an intriguing range of indirect evidence that this system has formed via accretion-induced collapse of a white dwarf, but the evidence at the present time is far from conclusive.

## ORCID iDs

Matteo Bachetti https://orcid.org/0000-0002-4576-9337
McKinley C. Brumback https://orcid.org/0000-0002-4024-6967

<footnote n="29">²⁹ Igoshev & Popov (2018) discussed timescales up to a few million yr.</footnote>
Felix Fürst https://orcid.org/0000-0003-0388-0560

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

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

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

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

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

References

Andersen, B. C., & Ransom, S. M. 2018, ApJL, 863, L13

Archibald, R. F., Kaspi, V. M., Tendulkar, S. P., & Scholz, P. 2016, ApJL, 829, L21

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., Harrison, F. A., Walton, D. J., et al. 2014, Natur, 514, 202

Basko, M. M., & Sunyaev, R. A. 1976, MNRAS, 175, 395

Begelman, M. C., King, A. R., & Pringle, J. E. 2006, MNRAS, 370, 399

Bildsten, L., & Brown, E. F. 1997, ApJ, 477, 897

Blaauw, A. 1961, BAN, 15, 265

Blackburn, J. K. 1995, in ASP Conf. Ser. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco, CA: ASP), 367

Blackburn, J. K., Shaw, R. A., Payne, H. E., Hayes, J. J. E. & HEASARC 1999, FTOOLS: A general package of software to manipulate FITS files, Astrophysics Source Code Library, ascl:9912.002

Boersma, J. 1961, BAN, 15, 291

Brightman, M., Harrison, F., Walton, D. J., et al. 2016, ApJ, 816, 60

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

Brightman, M., Walton, D. J., Xu, Y., et al. 2020, ApJ, 889, 71

Brumback, M. C., Hickox, R. C., Bachetti, M., et al. 2018, ApJL, 861, L7

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

Burderi, L., di Salvo, T., Menna, M. T., Riggio, A., & Papitto, A. 2006, ApJL, 653, L133

Canal, R., & Schatzman, E. 1976, A&A, 46, 229

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

Casella, P., Altamirano, D., Patruno, A., Wijnands, R. A. D., & van der Klis, M. 2008, ApJL, 674, L41

Conroy, C., Strader, J., van Dokkum, P., et al. 2018, ApJ, 864, 111

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

Dauser, T., Middleton, M., & Wilms, J. 2017, MNRAS, 466, 2236

Dib, R., & Kaspi, V. M. 2014, ApJ, 784, 37

Dib, R., Kaspi, V. M., & Gavriil, F. P. 2009, ApJ, 702, 614

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

Ergma, E. 1993, A&A, 273, L38

Fabrika, S., Ueda, Y., Vinokurov, A., Sholukhova, O., & Shidatsu, M. 2015, NatPh, 11, 551

Fragos, T., Linden, T., Kalogera, V., & Sklias, P. 2015, ApJL, 802, L5

Fragos, T., & McClintock, J. E. 2015, ApJ, 800, 17

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (3rd ed.; Cambridge: Cambridge Univ. Press)

Fryer, C., Benz, W., Herant, M., & Colgate, S. A. 1999, ApJ, 516, 892

Fuller, J. 2017, MNRAS, 472, 1538

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

Galloway, D. K., Morgan, E. H., Krauss, M. I., Kaaret, P., & Chakrabarty, D. 2006, ApJL, 654, L73

Galloway, D. K., Morgan, E. H., & Levine, A. M. 2004, ApJ, 613, 1164

Ghosh, P., & Lamb, F. K. 1979a, ApJ, 232, 259

Ghosh, P., & Lamb, F. K. 1979b, ApJ, 234, 296

Grisé, F., Kaaret, P., Pakull, M. W., & Motch, C. 2011, ApJ, 734, 23

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

HEASARC 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

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

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

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

Hut, P., & Paczynski, B. 1984, ApJ, 284, 675

Igoshev, A. P., & Popov, S. B. 2018, MNRAS, 473, 3204

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

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

Jiang, Y.-F., Cantiello, M., Bildsten, L., et al. 2018, Natur, 561, 498

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

Kalogera, V. 1998, ApJ, 493, 368

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

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

King, A. R. 2009, MNRAS, 393, L41

King, A. R., Pringle, J. E., & Wickramasinghe, D. T. 2001, MNRAS, 320, L45

Knigge, C., Coe, M. J., & Podsiadlowski, P. 2011, Natur, 479, 372

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

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

Kozai, Y. 1962, AJ, 67, 591

Kretschmar, P., Kreykenbohm, I., & Wilms, J. 1997, in Proc. 2nd INTEGRAL Workshop, The Transparent Universe, ed. C. Winkler, T. J.-L. Courvoisier, & Ph. Durouchoux (Noordwijk: ESA), 141

Lewin, W. H. G., van Paradigs, J., & van der Klis, M. 1988, SSRv, 46, 273

Lidov, M. L. 1962, P&SS, 9, 719

Luo, J., Ransom, S., Demorest, P., et al. 2019, PINT: High-precision pulsar timing analysis package, Astrophysics Source Code Library, ascl:1902.007

Maccarone, T. J., Kundu, A., Zepf, S. E., & Rhode, K. L. 2010, MNRAS, 409, L84

Manchester, R. N. 2018, in IAU Symp. 337, Pulsar Astrophysics the Next Fifty Years, ed. P. Weltevrede et al. (Cambridge: Cambridge Univ. Press), 197

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993

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

Middleton, M. J., Brightman, M., Pintore, F., et al. 2019a, MNRAS, 486, 2

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

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

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

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

Mohamed, S., & Podsiadlowski, P. 2007, in ASP Conf. Ser. 372, 15th European Workshop on White DWarfs, ed. R. Napiwotzki & M. R. Burleigh (San Francisco, CA: ASP), 397

Molkov, S., Lutovinov, A., Falanga, M., Tsygankov, S., & Bozzo, E. 2017, MNRAS, 464, 2039

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

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

Nomoto, K. 1987, ApJ, 322, 206

Pasham, D. R., & Strohmayer, T. E. 2013, ApJL, 774, L16

Perna, R., Bozzo, E., & Stella, L. 2006, ApJ, 639, 363

Pétri, J. 2016, JPIPh, 82, 635820502

Pfahl, E., Rappaport, S., Podsiadlowski, P., & Spruit, H. 2002, ApJ, 574, 364

Pike, S. N., Harrison, F. A., Bachetti, M., et al. 2019, ApJ, 875, 144

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

Pintore, F., Zampieri, L., Stella, L., et al. 2017, ApJ, 836, 113

Poelarends, A. J. T., Herwig, F., Langer, N., & Heger, A. 2008, ApJ, 675, 614

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

Ransom, S. 2011, PRESTO: Pulsar Exploration and Search TOLkit, Astrophysics Source Code Library, ascl:1107.017

Ransom, S. M. 2001, PhD thesis, Harvard Univ.

Ransom, S. M., Eikenberry, S. S., & Middleditch, J. 2002, AJ, 124, 1788

Rappaport, S., & Joss, P. C. 1977, Natur, 266, 683

Rappaport, S., Joss, P. C., & Webbink, R. F. 1982, ApJ, 254, 616

Ray, P. S., Guillot, S., Ho, W. C. G., et al. 2018, ApJ, 879, 130

Samus', N. N., Kazarovets, E. V., Durlevich, O. V., Kireeva, N. N., & Pastukhova, E. N. 2017, ARep, 61, 80

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Sci, 337, 444

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

Savitzky, A., & Golay, M. J. E. 1964, AnaCh, 36, 1627

Serim, M. M., Şahiner, Ş, Çerri-Serim, D., İnam, S. Ç., & Baykal, A. 2017, MNRAS, 471, 4982

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
Skinner, G. K., Bedford, D. K., Elsner, R. F., et al. 1982, *Natur*, **297**, 568

Taylor, J. H. 1992, *RSPTA*, **341**, 117

Taylor, J. H., & Weisberg, J. M. 1982, *ApJ*, **253**, 908

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

Usov, V. V. 1992, *Natur*, **357**, 472

Vasilopoulos, G., Haberl, F., Carpano, S., & Maitra, C. 2018, *A&A*, **620**, L12

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

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

Wang, Y.-M. 1995, *ApJL*, **449**, L153

Weisskopf, M. C., Brinkman, B., Canizares, C., et al. 2002, *PASP*, **114**, 1

Wiktorowicz, G., Lasota, J.-P., Middleton, M., & Belczynski, K. 2018, *ApJ*, **875**, 53