All at once: transient pulsations, spin down and a glitch from the Pulsating Ultraluminous X-ray Source M82 X-2

TL;DR

This paper reports the first detection of spin down and a glitch in the pulsating ultraluminous X-ray source M82 X-2, providing new insights into its timing behavior and challenging existing models.

Contribution

It presents the first observation of spin down and a glitch in a PULX, refining orbital parameters and expanding understanding of PULX timing properties.

Findings

Detected a spin down rate of 5x10^-11 Hz/s between 2014-2016

Identified the first positive glitch in a PULX system

Refined the orbital period of M82 X-2

Abstract

M82 X-2 is the first pulsating ultraluminous X-ray source (PULX)to be identified. Since the discovery in 2014, NuSTAR has observed the M82 field 15 times throughout 2015 and 2016. In this paper, we report the results of pulsation searches in all these datasets, 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 5x10^-11 Hz/s, which is in contrast to the strong spin up seen during the 2014 observations and represent 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 dataset. 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.

Tables

Table 1. Table 1: Pulsation search results from all available observations of M82 X-2.

ObsID	MJD	Date	LengthaaThe 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.	ExposureaaThe 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.	Count ratebbCount rates using FPMA + FPMB.	Meas. p.f.ccThe pulsed fraction is referred to the total X-ray flux of M82, since M82 X-1 and X-2 are not separable in NuSTAR.	Flux ratioddThe 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. in prep. when available, or from Brightman et al. (2019) otherwise (starred). Dashes indicate that a Chandra dataset was not available.	Corr. p.f.eeA 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.	$L_{39}$
			(ks)	(ks)	(s-1)	$\Delta F_{\mathrm{X2}}/F_{\mathrm{tot}}$(%)	$F_{\mathrm{X2}}/F_{\mathrm{X1}}$(%)	$\Delta F_{\mathrm{X2}}/F_{\mathrm{X2}}$(%)
80002092002	56681	2014-01-23	123	66	0.9	5	–	–	–
80002092004ffDetected only after MJD 56683.5.	56683	2014-01-25	171	90	1.0	4	–	–	–
80002092006	56686	2014-01-28	579	310	0.9	5	–	–	–
80002092007	56692	2014-02-04	562	306	0.9	7	132∗	12	8.1∗
80002092008	56699	2014-02-10	62	34	1.0	7	–	–	–
80002092009	56700	2014-02-11	213	115	0.9	9	–	–	–
80002092011	56720	2014-03-03	201	111	0.6	3	–	–	–
50002019002	57038	2015-01-15	56	31	0.7	6	40	27	5.7
50002019004	57042	2015-01-19	283	161	0.7	4	33∗	14	2.0∗
90101005002	57194	2015-06-20	56	37	2.3	3	32	17	19.8
80202020002	57414	2016-01-26	66	36	1.5	4	0∗	N.A.	0.0∗
80202020004	57442	2016-02-23	61	32	1.3	5	52	18	13.9
80202020006	57483	2016-04-05	54	31	0.8	6	44	26	6.7
30101045002	57493	2016-04-15	350	189	1.0	3	–	–	–
80202020008	57503	2016-04-24	67	40	1.1	5	79	13	10.8
30202022002	57543	2016-06-03	60	39	1.0	5	0∗	N.A.	0.1∗
30202022004	57571	2016-07-01	68	47	1.7	4	33	17	15.1
30202022008	57599	2016-07-29	67	43	1.6	4	4∗	N.A.	1.2∗
30202022010	57619	2016-08-19	69	44	1.2	4	16∗	37	4.3∗
90201037002	57641	2016-09-10	94	80	1.2	3	–	–	–
90202038002	57669	2016-10-07	71	45	0.9	5	5∗	N.A.	0.8∗
90202038004	57723	2016-11-30	68	43	0.8	5	0∗	N.A.	0.0∗
Det. puls.				1000 (58%)
Undet.				730 (42%)
Low flux				207

Table 2. Table 2: Updated Orbital parameters of M82 X-2, compared to the B14 fit and to a new fit of the 2014-only data. The quoted uncertainties represent approximate symmetrical 1- σ 𝜎 \sigma confidence intervals, as calculated by PINT . Values and error bars on a ​ sin ⁡ i 𝑎 𝑖 a\sin i and T 0 subscript 𝑇 0 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 40-d lever arm with respect to T 0 = subscript 𝑇 0 absent T_{0}= MJD 56682.0661, and including the local first derivative in the fit. The full 2014–2016 solution is based on the orbital phase measured in Fig.  3 and refined with PINT, fixing the first derivative to zero (see text).

Parameter	B14	2014-only fit	All data
$T_0$ (MJD)	56694.7327(1)	56682.0661(3)	56682.0661(3)
$a\sin i$ (lt-s)	22.225(4)	22.215(5)	22.215(5)
$P_{\mathrm{orb}}$ (d)	2.53260(5)	2.53287(5)	2.532948(4)

Table 3. Table 3: Spin parameters for the observations in 2014 using different models. Orbital parameters were fixed to the values in Table  2 . The reference epoch for all models is MJD 56683. 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 .

Parameter	Model 1	Model 2	Model 3	Model 4
F0 (Hz)	0.728566160(7)	0.72852171(3)	0.72855082(8)	0.72854408(11)
F1 (${10}^{-11}$ Hz/s)	0	4.758(3)	-2.98(2)	-5.29(3)
F2 (${10}^{-17}$ Hz/s2)	0	0	8.66(2)	10.93(3)
GL_F0 (Hz)	0	0	0	1.75(2)$\times {10}^{-5}$
GL_EP (MJD)	–	–	–	56685.8 (fixed)
r.m.s. ($\mu s$)	136	34	9.2	4.3
${\chi }^2$ ($\mu s^2$)	2698012	165709	12251	2704

Equations

$$L_{\mathrm{acc}}\approx L_{\mathrm{Edd}}\left(1+\ln \dot{m}\right),$$

$$L_{\mathrm{obs}}=L_{acc}/b\propto L_{\mathrm{Edd}}{\dot{m}}^2\left(1+\ln \dot{m}\right),$$

$$-\frac{\dot{p}}{p}\propto pn({\omega }_{\mathrm{fast}})$$

$$n({\omega }_{\mathrm{fast}})\approx \frac{7/6-(4/3){\omega }_{\mathrm{fast}}+(1/9){\omega }_{\mathrm{fast}}^2}{1-{\omega }_{\mathrm{fast}}}.$$

$$-\frac{\dot{p}}{p}\propto p$$

$$-\frac{\dot{p}}{p^2}=\dot{\nu }\propto L_{37}^{6/7}n({\omega }_{\mathrm{fast}}),$$

$$-\frac{\dot{p}}{p^2}=\dot{\nu }\propto L_{37}^{3/7}n({\omega }_{\mathrm{fast}}).$$

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

$$\begin{gathered} \frac{\dot{a}}{a}=\frac{2{\dot{P}}_{\mathrm{orb}}}{3P_{\mathrm{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] \end{gathered}$$

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

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

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

$${\dot{P}}_{\mathrm{orb}}$$

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

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

Full text

stars: neutron — pulsars: individual (M82 X-2) — X-rays: binaries — ULX

††facilities: NuSTAR (Harrison et al., 2013), Chandra (Weisskopf et al., 2002), ATNF pulsar catalogue (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 (Price-Whelan et al., 2018), PINT (Luo et al., 2019), PRESTO (Ransom, 2011)

®

1 Introduction

Ultraluminous X-ray sources (ULX) 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 ultraluminous X-ray sources (PULX111also referred to as ultraluminous x-ray pulsars, ULXP) 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 average222For a detailed view see Pintore et al. 2017; Walton et al. 2018; Koliopanos et al. 2017. 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 stars333For deeper investigations of ULX populations, see Middleton & King (2017); Wiktorowicz et al. (2018).

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 observations 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; Perna et al., 2006a; Burderi et al., 2006).

Accretion onto magnetized neutron stars happens through the interaction of the accretion disk with their strong magnetic field. The general framework of accretion onto magnetized objects is described in a number of papers (e.g. Ghosh & Lamb, 1979a, b; Wang, 1995) and reviewed in Frank et al. (2002). The disk is truncated at a distance $R_{\rm in}$ (often called the Alfvén radius) from the NS, where magnetic stresses equal the “ram pressure” from the matter in the disk. Therefore, the same $R_{\rm 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_{\rm in}$, also penetrates the disk beyond the inner radius. The interaction between the matter and the magnetic field at a given distance from the NS produces a torque on the star that depends on the relative angular velocity of the matter with respect to the star. Calling $R_{\rm 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_{\rm 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 in this theoretical framework with good success. As a final ingredient in standard accreting pulsars, the X-ray luminosity is considered a proxy of 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 spherization radius $R_{\rm sph}$. In super-Eddington disks, inside $R_{\rm 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

[TABLE]

where $\dot{m}=\dot{M}/\dot{M}_{\rm Edd}$. However, this is probably not the bolometric luminosity measured by the observer, because radiation might be beamed. King (2009) describe 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

[TABLE]

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

If $R_{\rm sph}$ >$$R_{\rm in}, the disk becomes locally super-Eddington before reaching the inner radius. So, the actual mass accretion rate (the mass that reaches $R_{\rm in}$ and accretes on the NS) 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., 2016). 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_{\rm 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_{\rm 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 onto 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, constrain the total mass exchange in the system. The observation was performed on UT 2016-04-15 – 2016-04-19 (MJDs 57493.29 – 57497.34).

NuSTAR does not resolve M82 X-2 from another ULX, M82 X-1, only 5$\arcsec$ 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, Brightman et al. in prep.).

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, 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 development444The temperature-driven clock correction is now used to produce the official clock correction files with NuSTAR. The new clock correction is default since barycorr v2.2. Bachetti, Markwardt et al. in prep, and used PINT555www.github.com/nanograv/PINT for the satellite orbit calculations. The pulse period of M82 X-2 is three 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 nupipeline666During 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 account properly for leap seconds.. 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 and not just from a single source; we removed the relevant intervals by modifying the good time intervals (GTI). Following B14, we selected photons from a 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_{\rm x,max}\sim 10^{41}$ erg/s) and M82 X-2 ($L_{\rm x,max}\sim 2\times 10^{40}$ erg/s) (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) and Brightman et al. sub.

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 HENDRICS777www.github.com/StingraySoftware/HENDRICS (Bachetti, 2018), Stingray 888www.github.com/StingraySoftware/stingray (Huppenkothen et al., 2019, 2016) PINT (Luo et al., 2019), and PRESTO (Ransom, 2011, 2001).

In the 2-Ms long campaign performed in 2014, where the pulsar was first detected, pulsations did not have a constant r.m.s. A large variation of pulsed fraction was reported by B14 and is also shown in Fig. 9. 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 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 two 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.

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 HENDRICS999 Rather than calculating the $Z^{2}_{n}$ statistics on the single events, this software pre-folds 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).

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 (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 8 times finer than the standard (e.g. 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 on “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 $\sim$200 s in order to account for a possibly imperfect orbital solution101010200 s was chosen as it is approximately the error on $T_{0}$ that would produce a shift of $\sim$half a pulsar rotation over an orbit. 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 in segments of 30 ks, with sliding windows overlapping by a factor $\sim 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_{\rm orb}$). This was clearly seen in the 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_{\rm orb}$.

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

We first used the accelsearch tool, that performs an accelerated search of pulsations based on the FFT and matched filters. We used the lowest detection limit (-sigma=1 on the command line), and specified that 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 beats of fundamental frequencies of the NuSTAR detector111111e.g. a small dead time component that repeats with a frequency of $\sim 890$ Hz in datasets taken in CP mode, usually visible in very bright sources or to 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). prepfold 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 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 ten 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 if a good candidate was consistently found with similar starting periods) and 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 the 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 a $\sim$1.37 s pulsation, Doppler shifted by orbital motion with $P_{\rm orb}\sim$2.52 d and $a\sin i/c\sim 22$ lt-sec, the orbital parameters known from B14. We refined the spin solution using the full observation (instead of a 80-ks segment) using HENphaseogram121212The signal-to-noise in single intervals was not adequate to make more sophisticated timing using, e.g. PINT or tempo2 at this step. These local spin derivatives give effects two orders of magnitude above the reasonable interval for the spin-up of the source, which is $|\dot{\nu}|\lesssim 2\times 10^{-10}$ Hz/s. 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 (Fig. 3). We can expect intrinsic (torque-driven) spin-up or spin-down whose magnitude is up to some $10^{-10}$ s/s. The spin-up parameters and the orbital parameters are degenerate, and we include these considerations when calculating the improved error bar on the orbital period $P_{\rm orb}$.

Nonetheless, the 2.5-year lever arm between this new measurement and the original 2014 dataset allows for the error bar on $P_{\rm orb}$ to be significantly reduced (Table 2). With this rough orbital solution, we were able to fold the events in sub-intervals of the new observation and calculate the arrival times of the pulses with the fftfit method (Taylor, 1992) as implemented in HENphaseogram. Then, we used PINT to refine the orbital and spin parameters, fitting for the frequency, the first derivative and the orbital period. Given the small observation length, the spin derivative is not constrained by the fit, therefore we fixed it at 0. Finally, we found a well-constrained solution and an improved value of $P_{\rm 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_{\rm 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 $\sim$0.0001 d, 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 Fig. 5, 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 (Fig. 9) of $\dot{\nu}\sim-5.8\times 10^{-11}$ Hz/s 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^{2}_{1}$ (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 non-detections. To do so, for each dataset, we maintained the GTIs and the total number of photons in each interval and randomized the time of the events inside each GTI. There is no significant broad-band 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 sub-intervals of the observations, and aligning the folded sub-profiles differently for slightly different values of the frequency (similarly 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 sub-intervals 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 1) that the number of bins in the profile $N$ is large enough that the pulse shape features are distinct in different bins and 2) that the number of profiles that are being shifted is at least twice the maximum shift of the sub-profiles. We verified that calculating these statistical tests this way does not depart significantly from the expected statistics (More details in the Appendix and Fig. 12). Note that this new folding algorithm, as of August 13th, 2019, is 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 (Fig. 3). We only find a hint of pulsations in ObsID 30101045002, close to the predicted spin-down line and to the detection limit, that we plot in Fig. 9 with a grey 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 why 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 ObsIDs 800020920XX 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.

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 ($\ddot{\nu}$) 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 plugin of TEMPO2, we calculate two local spin derivatives (F0 and F1) from the pulse times of arrival. Fig. 7 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{\nu}$ necessary to produce such a large frequency increase on timescales of  1 day is $\sim 10^{-10}$ Hz/s. Based on simple arguments from accretion theory, such a strong frequency derivative would need a significant change of mass accretion rate and, presumably, of 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 (Bachetti et al. in prep.), 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 pulsed frequency, but we verified that the scatter of the offset measurements from the Malindi station is $\leq$100 $\mu$sec 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 non-detection 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 anti-glitches 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): the authors find that, despite a large drop of luminosity of the source NGC 300 ULX1, 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 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 $\sim 60$ days. For the physical properties of this system, however, this is problematic. While it would be easy to produce a a periodically-changing mass transfer rate if the period corresponded to the orbital period, this is far from obvious for a super-orbital periodicity. Assuming Roche Lobe overflow, possible mechanisms to increase periodically the pressure at the L1 Lagrangian point might involve additional orbiting objects, like the so-called Kozai mechanism (Lidov, 1962; Kozai, 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 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 $\beta$-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-d 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) rule 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 non-detections

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 r.m.s. $r$, and only weakly by the observing time $T$. Lewin et al. (1988) show that the significance of pulsations goes roughly as $I^{2}/(I+B)r^{2}T^{1/2}$. Looking into detail at the table, we see that 2014 observations were executed in a very favorable situation, with M82 X-2 brighter than M82 X-1 during ObsID 80002092007, a high (even if not at the highest values) 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.

Fig. 10 shows that the increase of 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 to intepret the increasing torque implied by the strong second derivative $\ddot{\nu}$ 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, 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 two weeks before and after, 30202022010 and 90202038002, despite the overall NuSTAR count rate being similar131313Note that change of flux from M82 X-2 on these time scales is compatible with the super-orbital period of 60 days shown by Brightman et al. (2019)..

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 is 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. M82 X-2 is not the only ULX showing intermittent pulsations. Recently, Sathyaprakash et al. (2019) found intermittent pulsations in NGC 1313 X-2. ULXs are peculiar objects, probably Roche-lobe filling like low-mass X-ray binaries (that have much lower companion masses) but with companion masses more similar to high-mass X-ray binaries (that mostly accrete via wind)141414For completeness, we mention that Mellah et al. (2019) proposed an intermediate Roche-wind accretion model, initially proposed by Mohamed & Podsiadlowski (2007), for ULXs.

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^{9}$G) and sub-Eddington. Casella et al. (2008) reported pulsations lasting only $\sim 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). A0538-66, an HMXB, showed pulsations in 1982 (Skinner et al., 1982) and then never since (see Kretschmar et al., 2004, 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 (r.m.s., 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 non-pulsed and the 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. Pinto et al., 2016; Fabrika et al., 2015; 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) model long super-orbital periods observed in ULXs as the precession from a disk launching a conical wind. Precession is observed in an archetipal super-critical source, SS433, 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. in prep. do not find a 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 pile-up, 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 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_{\rm 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 a relation

[TABLE]

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

[TABLE]

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

[TABLE]

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 (Fig. 11). 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 Ghosh & Lamb, 1979b)

[TABLE]

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

[TABLE]

Again, Fig. 11 shows clearly for both models that the spin behavior of M82 X-2 in 2014 is inconsistent with $\omega_{\rm fast}<<1$, while it can in principle be reproduced if we assume $\omega_{\rm fast}\sim 1$. Therefore, the basic assumption of most concurring models on M82 X-2 (e.g. Bachetti et al., 2014; Dall’Osso et al., 2015; Tsygankov et al., 2016; King et al., 2017) Here we calculate spin derivatives from a robust interpolation151515We use the method by Savitzky & Golay (1964), as implemented in scipy.signal 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 reported for LMC X-4 by Molkov et al. (2017), where an alternating spin-up and down behavior was observed over many decades. In that case, one of the possible explanations was a Recycling magnetosphere model (Perna et al., 2006b) where accretion is only possible at certain spin phases.

Independently from accretion torque, the pulsar might be spinning down because of its own dipole radiation. According to the classic magnetic dipole radiation formula161616See any textbook on pulsars, e.g. https://www.cv.nrao.edu/course/astr534/Pulsars.html, if the spin-down of the pulsar is given by dipole radiation, we can estimate the magnetic field as:

[TABLE]

where $I$ is the moment of inertia of the neutron star and $R$ its radius. Let us assume a fine tuned equilibrium from the accretion torque, with no net spin-up or down from accretion and the spin-down dominated by dipole emission. Substituting the spin period and the 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. This would be in line with the estimate made by Tsygankov et al. (2015), interpreting the high- and low-luminosity states of the pulsar as alternating accretion-propeller regimes and assuming that there is no mass loss before accretion onto the pulsar. Future robust measurements of the total mass exchange will give a better estimate of the position of the spherization radius, and will help to test the predictions of these and other competing models.

4.4 Glitch: strong spin-down before 2014?

There are very few glitches reported from accreting neutron stars. KS 1947+300 showed a large glitch of $\Delta\nu/\nu\approx 3.5\times 10^{-5}$ in RXTE observations in the early 2000s (Galloway et al., 2004), comparable with $\Delta\nu/\nu\approx 10^{-5}$ reported in Section 3.5. Similarly, Serim et al. (2017) report on a strong glitch from the slow accreting pulsar SXP1062, 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, and notably the ultraluminous pulsar NGC 300 ULX-1 (Ray et al., 2018), show *anti-*glitches instead – sudden slow-downs of the pulsar. A non-magnetar glitch is believed to be due to the independent rotation of the NS core with respect to the rest of the NS. A neutron star contains superfluid material in the core and in the 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 the 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 anti-glitch 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 in the anti-glitches likely represent the effect of reconnection between the faster crust and slower core. It is possible that the glitch we observe at epoch 56685.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 NS 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 dataset. 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 three years 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}_{\rm orb}\sim-10^{-7}$ d/d implying an orbital evolution time-scale of 30000 yrs and a mass exchange almost three orders of magnitude above the Eddington limit (see below).

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

[TABLE]

where $M_{c}$ an $\dot{M}_{c}$ are the mass and mass loss rate from the companion (negative), $M_{p}$ an $\dot{M}_{p}$ the same quantities for the pulsar, $q=M_{p}/M_{c}$, $\beta=-\dot{M}_{p}/\dot{M}_{c}$, and $\alpha$ the specific angular momentum of the lost mass in units of $2\pi a^{2}/P_{\rm orb}$. In the conservative scenario, $\beta=1$ and (9) becomes

[TABLE]

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

[TABLE]

and the (9) reduces to

[TABLE]

which is a relatively small correction to (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

[TABLE]

This $\dot{P}_{\rm orb}$ corresponds to a change of orbital period of $\sim 1$ s/yr, and produces a shift of the orbital phase by a few 100 s over a few years. Again, since $\dot{M}_{c}$ is negative, also this $\dot{P}_{\rm orb}$ is 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, the remaining extreme scenario where most mass is lost from the system in form of winds from the donor star, we have

[TABLE]

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

It is likely that in a real-life scenario both phenomena should happen, stabilizing the mass transfer over long time scales. Fragos et al. (2015) discuss 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 two years, 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 that:

• •

The pulsar alternates phases with a large spin-up to 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.

• •

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.

• •

The neutron star has very strong glitches, probably due to the rapid spin evolution. The observed positive glitch suggests a strong spin-down to have 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.

MB thanks the Fulbright Visiting Scholar Program for supporting a nine-month visit at Caltech. MM appreciates support from an STFC Ernest Rutherford Fellowship. Part of the pulsar search and optimization software was developed in the framework of CICLOPS – Citizen Computing Pulsar Search, a project supported by POR FESR Sardegna 2014 – 2020 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, Wynn Ho for pointing out a miscitation 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.

Appendix A On the detection limits used in the pulsation search.

Fig. 12 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 data, and on (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 follows very closely 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 alter significantly 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 from the analysis technique used.

To evaluate the $p$-value, in order not to increase synthetically 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, 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):

PSR M82_X2_OBSID_007_009 RAJ 9:55:42.14400000 DECJ 69:40:26.00400000 F0 0.7285667095459839939 1 8.741567555166416003e-08 F1 6.4250081648893346535e-11 1 4.9287115294950599775e-13 F2 2.2525682844363675388e-16 1 7.8038412567833659266e-18 F3 -5.939687798010528739e-22 1 5.0019855875541822192e-23 PEPOCH 56695.507895798600000 PLANET_SHAPIRO N BINARY BT PB 2.5328840644740456124 1 0.00016109793561046761965 A1 22.204213222811447 1 0.012032723626796296 ECC 0.0 T0 56682.068219262420266 1 0.001018781929672351325 OM 0.0

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_{\rm orb}$:

PSR M82_X2_OBSID_011 RAJ 9:55:42.14400000 DECJ 69:40:26.00400000 F0 0.72875923701101958546 1 3.364666714127084281e-07 F1 1.189117392632877145e-10 1 6.9989856290673740436e-12 PEPOCH 56720.878600229010000 PLANET_SHAPIRO N BINARY BT PB 2.5328387434901325451 1 4.4959856101332736067e-05 A1 22.204213222811447 0 0.012032723626796296 ECC 0.0 T0 56682.068219262420266 0 0.001018781929672351325 OM 0.0

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:

PSR M82_X2_OBSID_90201037 RAJ 9:55:50.20800000 DECJ 69:40:46.99200000 F0 0.72390599836134799007 1 7.123714905580943099e-07 PEPOCH 57641.998517226645995 PLANET_SHAPIRO N BINARY BT PB 2.5329477619141624702 1 4.0457956366063808576e-06 A1 22.204213222811447 0 0.012032723626796296 ECC 0.0 T0 56682.068219262420266 0 0.001018781929672351325 OM 0.0

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 magnetar-scale 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 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 10000 years171717Igoshev & Popov (2018) discuss timescales up to a few million years, but that 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, which 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 accretion inducted 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 collapse 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 the glitch (and as previously suggested for this source in the literature, e.g. Tsygankov et al. 2015, Mushtukov et al. 2015, Ekşi et al. 2015, Dall’Osso et al. 2015), 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.