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

# The Slowest Spinning X-Ray Pulsar in an Extragalactic Globular Cluster

Ivan Yu. Zolotukhin<sup>1,2,3</sup>, Matteo Bachetti<sup>4</sup>, Nicola Sartore<sup>5</sup>, Igor V. Chilingarian<sup>2,6</sup>, and Natalie A. Webb<sup>1,5</sup>

Received 2016 February 15; revised 2017 March 3; accepted 2017 March 20; published 2017 April 25

## Abstract

Neutron stars are thought to be born rapidly rotating and then exhibit a phase of rotation-powered pulsations as they slow down to 1–10 s periods. The significant population of millisecond pulsars observed in our Galaxy is explained by the *recycling* concept: during an epoch of accretion from a donor star in a binary system, the neutron star is spun up to millisecond periods. However, only a few pulsars are observed during this recycling process, with relatively high rotational frequencies. Here we report the detection of an X-ray pulsar with $P_{\text{spin}} = 1.20$ s in the globular cluster B091D in the Andromeda galaxy, the slowest pulsar ever found in a globular cluster. This bright (up to 30% of the Eddington luminosity) spinning-up pulsar, persistent over the 12 years of observations, must have started accreting less than 1 Myr ago and has not yet had time to accelerate to hundreds of Hertz. The neutron star in this unique wide binary with an orbital period $P_{\text{orb}} = 30.5$ hr in a 12 Gyr old, metal-rich star cluster accretes from a low-mass, slightly evolved post-main-sequence companion. We argue that we are witnessing a binary formed at a relatively recent epoch by getting a ~0.8 $M_{\odot}$ star in a dynamical interaction—a viable scenario in a massive, dense globular cluster like B091D with high global and specific stellar encounter rates. This intensively accreting non-recycled X-ray pulsar therefore provides a long-sought missing piece in the standard pulsar recycling picture.

**Key words:** astronomical databases: miscellaneous – galaxies: individual (M31) – globular clusters: individual (Bol D91) – pulsars: individual (3XMM J004301.4+413017) – virtual observatory tools – X-rays: binaries

## 1. Introduction

Around 2000 pulsars are known, where the majority of these are “regular” pulsars that have pulses periods between tens of milliseconds to approximately a second, and magnetic field strengths of ~10<sup>12</sup> G (Manchester et al. 2005). These pulsars show a general spin down to longer periods. However, a few hundred of these pulsars show much shorter periods, of the order of a millisecond, along with lower magnetic fields of ~10<sup>8</sup> G. It is believed that these millisecond pulsars (MSPs) are the descendants of neutron stars/pulsars found in X-ray binaries. Accretion onto a neutron star from a close companion is believed to transfer angular momentum to the neutron star, spinning it up to periods of milliseconds (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). This is strongly supported by both the discovery of an MSP in an X-ray binary system (SAX 1808.4-3658, Wijnands & van der Klis 1998) and the presence of kilo-Hertz Quasi-Periodic Oscillations in many LMXBs, which have been found to show millisecond pulsation periods (see van der Klis 1998 and references therein). More recently, Archibald et al. (2009) showed that the previously accreting millisecond pulsar FIRST J102347.67+003841.2 had ceased to accrete and radio pulsations could subsequently be observed, thus supporting the “recycled” pulsar idea (see Patruno & Watts 2012 for a review).

Many of the known MSPs are found in Galactic globular clusters (GCs), where as of the end of 2015, almost 140 MSPs have been detected.<fnref n="7" /> GCs are dense spherical systems of ~10<sup>4</sup>–10<sup>6</sup> old stars (e.g., Hénon 1961). Their old age implies that they

should also contain many compact objects (e.g., Hut et al. 1992). Stellar encounters, which are extremely rare in lower density regions, can occur in GCs on timescales comparable with or less than the age of the universe. This would indicate that many Galactic GC stars have undergone at least one encounter in its lifetime. Encounters between stars is one way in which binaries can be produced. The encounter rate (Γ) due to tidal capture (Fabian et al. 1975) is proportional to the encounter cross-section, the relative velocity of the stars, and the number density of stars in the cluster (core). Both primordial binary systems and those formed due to encounters should exist in GCs, due to the dense environments, but encounters between a binary and either a single star or a binary system would more readily occur as the cross-sections are significantly larger, thus increasing the likelihood of an encounter. This explains the large number of recycled pulsars that we observe in GCs.

However, a small number (6) have periods greater than 0.1 s, and these pulsars have not yet been fully recycled. The longest of these, B1718-19 in NGC 6342, has a period of 1.004 s and a magnetic field of ~10<sup>12</sup> G, typical of a “regular” pulsar. It appears to be a young pulsar, with a characteristic age of 1 × 10<sup>7</sup> years (Lyne et al. 1993). Lyne et al. (1993) propose that this young pulsar originated either from a collision between an old neutron star and a cluster star, or a white dwarf accreting from a companion underwent an accretion-induced collapse (Michel 1987). Verbunt & Freire (2014) suggest that the current main-sequence companion of an old neutron star in that system has replaced the original one in an exchange encounter.

<footnote n="1"><sup>1</sup>Université de Toulouse; UPS-OMP, IRAP, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="2"><sup>2</sup>Sternberg Astronomical Institute, Moscow State University, Universitetskij pr., 13, 119992, Moscow, Russia</footnote>

<footnote n="3"><sup>3</sup>Special Astrophysical Observatory of the Russian Academy of Sciences, Nizhniy Arkhyz 369167, Russia</footnote>

<footnote n="4"><sup>4</sup>INAF/Osservatorio astronomico di Cagliari, via della Scienza 5, I-09047 Selargius, Italy</footnote>

<footnote n="5"><sup>5</sup>CNRS, IRAP, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="6"><sup>6</sup>Smithsonian Astrophysical Observatory, 60 Garden St. MS09, Cambridge, MA 02138, USA</footnote>

<footnote n="7">⁷ http://www.naic.edu/~pfreire/GCpsr.html</footnote>
Following the recent detection of coherent pulsations from
an ultraluminous X-ray source (ULX) in the M82 galaxy
(≈3.5 Mpc away) using NuSTAR data (Bachetti et al. 2014),
which showed that this bright source was in fact a neutron star,
we started to search archive XMM-Newton data to find similar
sources, in order to address questions such as how can such
super-Eddington luminosities be possible in a neutron star (see,
e.g., Mushtukov et al. 2015; King & Lasota 2016, for further
discussion).

In this paper, we describe the analysis that made it possible
to detect a 1.2 s pulsar in the X-ray binary 3XMM J004301.4
+413017, associated with the GC B091D from the Revised
Bologna Catalog of M31 GCs (RBC V.5; Galleti et al. 2004),
using public data. For brevity we denote the pulsar XB091D,
after its host GC designation. This is the first persistently
accreting X-ray pulsar hosting a neutron star detected in M31,<fnref n="8" />
and it also has the longest period among all known pulsars
(rotation-powered and accreting) in GCs, being more than an
order of magnitude slower than the mildly recycled accreting
pulsar from the Terzan 5 GC (Papitto et al. 2011). We note that
this pulsar has also recently been detected by Esposito et al.
(2016), but these authors interpret the nature of this source
quite differently. In this paper, we discuss possible evolu-
tionary scenarios that may produce such slowly rotating
neutron stars in GCs.

The content is organized as follows: Section 2 briefly
describes the *XMM-Newton* photon database (which will be
described in detail elsewhere) and the data set used for the
initial pulsation detection, as well as the *XMM-Newton* pulsar
factory analysis methods which resulted in the automated
detection of the pulsed X-ray emission; Section 3 covers the
manual blind search for pulsed emission in all available *XMM-
Newton* data for this source and the determination of the
source’s orbital parameters; Sections 4 and 5 summarize our
findings on the timing and spectral properties of this source,
which are then discussed in Section 6 where we also argue on
the possible origin and evolution of this system.

## 2. Pulsations Search Data and Methods

### 2.1. Photon Database

For this study, we created a database of all photons registered
by the EPIC/pn detector of the *XMM-Newton* satellite while
operating in imaging mode during the 7781 observations that
took place between 2000 and mid 2013. These are the same
observations that were used to compile the 3XMM-DR5 catalog
(Rosen et al. 2015). The photon data were taken from the event
lists called *PIEVLI* files publicly available from the *XMM-
Newton* science archive.<fnref n="9" /> These are science-ready data products
that come from the pipeline run by the *XMM-Newton* Survey
Science Centre (XMM SSC).<fnref n="10" /> Each file is a binary table with
one row per event that contains the following information:
(a) event time (satellite clock), (b) raw CCD pixel and projected
sky coordinates relative to the nominal pointing position,
(c) corrected and uncorrected event energy, (d) event quality

flag, and (e) CCD number. Each PIEVLI file also contains good
time intervals determined separately for each pn CCD by the
XMM SSC pipeline software. We point out that it is not
necessary to retrieve the complete XMM-Newton observation
archive (ODF archive) for the massive scale timing analysis with
XMM-Newton.

These event lists represent the lowest and most detailed data
level we use for initial timing studies. These are accompanied
by higher level data in the form of the 3XMM-DR5 catalog of
X-ray sources and their individual observations, acting
essentially as an index for the navigation within the photon
database. The connection between these two very different sets
of information is achieved by means of the known transforma-
tion between pixel coordinates of the events and world
coordinates of the X-ray sources from the 3XMM-DR5 catalog.
We used the World Coordinate System (WCS) transformation
data available for every EPIC/pn exposure in the PIEVLI file.
In this fashion, one can make an event extraction for any
detection from the 3XMM-DR5 catalog using its world
coordinates, right ascension, and declination.

The catalog of X-ray sources was cross-matched with other
catalogs of astrophysical objects in order to find possible
counterparts and determine the source type or location within a
galaxy. This gives us the possibility to easily extract arbitrary
photon lists of, e.g., all photons coming from the M31 galaxy
or all photons from known magnetars.

The last stage before being able to launch our timing analysis
codes over the sequence of extracted photon lists is the
barycentric correction, i.e., correcting the event times to the
solar system barycenter. For this purpose, we used the *XMM-
Newton* orbit files fed to the *barycen* task from the *XMM-
Newton* SAS version 13.5 software and the source positions
from the 3XMM-DR5 catalog.

We note that it is now possible to access this photon database
as well as the 3XMM-DR5 catalog data through a convenient
Web interface<fnref n="11" /> that we developed while working on the
3XMM-DR5 catalog compilation and this project. In particular,
one can extract barycentered photons in arbitrary regions from
the observation level event lists using nothing but a Web
browser. More details on this Website for the quick-look
science analysis of the *XMM-Newton* data will be presented in a
separate paper (I. Zolotukhin et al. 2017, in preparation).

### 2.2. Pulsar Factory Analysis Software

The photon database provides a way to easily extract
calibrated and barycentered event lists for arbitrary sets of
astrophysical sources observed with the EPIC/pn camera on
board the *XMM-Newton* satellite. We developed an analysis
software aimed at finding coherent pulsations in the photon
database, optimized for high throughput.

It consists of few analysis layers. First, for each event list, it
produces the power density spectrum (PDS) adapting the time
binning so that the Nyquist frequency is well above 2 kHz and
the total number of bins between gaps in the light curve is
factorable with small prime numbers for better performance
with the * numpy * FFT algorithm. In the PDS, the algorithm
searches for peaks exceeding the 99% detection level,
including the number of trials, following the standard rules
from Leahy et al. (1983). If such peaks are found, it launches
the * Zn*² test (Buccheri et al. 1983) in the vicinity of their

<footnote n="8">⁸ Other known X-ray pulsating sources in M31 are a transient pulsar candidate XMMU J004415.8+413057 in a high-mass X-ray binary with a period of 197 s (Trudolyubov et al. 2005), and two supersoft sources powered by accreting white dwarfs: XMMU J004252.5+411540 with a period of 217.7 s (Trudolyubov & Friedhorsky 2008) and XMMU J004319.4+411759 with a period of 865.5 s (Osborne et al. 2001).</footnote>

<footnote n="9">⁹ http://xmm.esac.esa.int/xsa/</footnote>

<footnote n="10">¹⁰ http://xmmssc.irap.omp.eu</footnote>

<footnote n="11">¹¹ http://xmm-catalog.irap.omp.eu</footnote>
frequencies. If a peak from the PDS is confirmed with the $Z_n^2$ test, an accelerated search with the PRESTO code (Ransom 2001) is launched to further constrain or discard this pulsation frequency candidate. As the final step of the automated procedure, the pulsar factory checks this frequency in other observations of this source from the 3XMM-DR5 catalog. If the detection is highly significant (>nσ, where n can vary between runs) or the same frequency is found in other detections, an operator is notified.

We first tested the XMM-Newton pulsar factory code using XMM-Newton data on a set of known X-ray-bright magnetars from the McGill Online Magnetar Catalog<fnref n="12" /> (Olausen & Kaspi 2014). After achieving stable detection of known coherent pulsation periods in this automated regime, we launched a larger period search on the unstudied sample of X-ray sources in the catalog. The M31 galaxy was an obvious choice of survey region for its large number of known X-ray sources and its relatively low distance, and was included in our first survey.

### 3. Detection and Orbital Parameters

During the M31 test run, the automated XMM-Newton pulsar factory algorithm detected pulsations, at about the same period, in three observations (ObsIDs: 0112570101, 0505720301, P0650560301) of the source 3XMM J004301.4+413017.<fnref n="13" /> The $Z_n^2$ test triggered by the detection confirmed the candidate. (see Figure 1) These detections are fully reproducible online from the source’s event extraction pages, e.g., for observation 0112570101: http://xmm-catalog.irap.omp.eu/pievli/101125701010068.

In order to look for the detected pulsation in more ObsIDs, we ran an accelerated search with PRESTO in all ObsIDs containing the source. The very high values of the period derivative required by the accelerated search, the clear improvement of detection significance when adding a second derivative in the search, and the shape of the track in the phaseogram shown by PRESTO, pointed to the presence of orbital modulation. Largely following the same procedure described in Bachetti et al. (2014), we cut the two longest ObsIDs into chunks, 10–30 ks long, and ran an accelerated search with custom-made software, searching the solution in the $\nu-\dot{\nu}$ plane (where $\nu$ indicates the pulse frequency) that yielded the highest $Z_n^2$ statistics.

Inside the long ObsIDs (0112570101, 0690600401), the best-solution frequency and frequency derivative clearly followed a sinusoidal law with a period between 1 and 2 days, as expected from orbital modulation. We fitted simultaneously the two values of frequency and frequency derivative in ObsID 0690600401 with the expected variation due to orbital motion, and obtained a first estimate of the orbital parameters ($P_{\text{orb}} \sim 1.2$ days, $a \sin i/c \sim 2.60$ 1s, $T_{\text{asc}} \approx \text{MJD}56104.789$, where $a \sin i/c$ is the projected semimajor axis and $T_{\text{asc}}$ the time of passage through the ascending node). Starting from the first rough estimate of these parameters and by trial and error trying to align the pulses in the phaseogram in the chunks first, then calculating TOAs with a custom implementation of the `fftfit` method (Taylor 1992) and using TEMPO2 to fit an orbital solution with the ELL1 model,<fnref n="14" />

<table><caption>**Table 1**
Best Orbital Solution Found in this Paper</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>P<sub>orb</sub></td><td>1.27101304(16) days</td></tr><tr><td>T<sub>asc</sub></td><td>MJD 56104.791(26)</td></tr><tr><td>a sin i/c</td><td>2.89(13)1s</td></tr><tr><td>e</td><td>&lt;0.003</td></tr></tbody></table>

Note. Brackets indicate $1\sigma$ uncertainties, as returned by TEMPO2. The upper limit on the eccentricity is the maximum uncertainty returned by the ELL1 model, for values of the eccentricity always consistent with 0.

we reached a solution valid for ObsID 0690600401. We then applied the solution to the other long ObsID, refining the orbital parameters so that every residual orbital modulation was eliminated by assuming a constant spin through the observation. The solution found in this way is the following: $P_{\text{orb}} \approx 1.2695$ days, $a \sin i/c \approx 2.886$, and $T_{\text{asc}} \approx 56104.7907$. Then, we addressed the ObsIDs in between, using “quantized” values of the orbital period that conserved the ascending node passages close to the two long ObsIDs. For every value, we calculated the scatter that it produced on the TOAs and selected the one that produced the lowest scatter. The eccentricity fitted by ELL1 was always consistent with 0 at the ~2σ level, with an uncertainty of 0.002–0.003. We use this last number as an upper limit on the eccentricity. The full solution is in Table 1.

We used this solution to look for pulsations in the remaining ObsIDs and refine the estimate. In the next section we describe this procedure.

### 4. Refined Timing Analysis

Except for a few cases (ObsIDs: 0112570101, 0650560301, 0690600401), all other individual observations yielded photon statistics that were too poor to determine the pulsation period and the pulse shape with enough statistical significance for detailed interpretation.

Therefore, we attempted the search for coherent pulsations by combining several data sets spanning two to five months in different years of observations by assuming that the period did not change among individual observations within each block. First, we corrected all the photon arrival times using the orbital elements of the binary system reported above. Then we used the period grid search using the S statistics with regularization (see Appendix for details), leaving the orbital phase as an additional free parameter.

Then, for each year, we started with the first observation typically taken in late December or early January. We then started adding observations, checking that the S statistics around the probable period, according to the increase in the exposure time, suggest that the pulsations are still coherent. If the S statistics (Leahy et al. 1983) did not increase or decrease when we added the additional observation, we concluded that the period changed significantly and started a new block. The results of our timing analysis are provided in Table 2, and the blocks of observations used for coherent searches are separated by horizontal lines. We estimated uncertainties of the period measurements analytically from the photon statistics, exposure time, and the pulse properties as explained in the Appendix. In Figure 2, we provide the nine recovered pulse profiles between 2002 and 2012.

<footnote n="12">¹² http://www.physics.mcgill.ca/~pulsar/magnetar/main.html</footnote>

<footnote n="13">¹³ See the source Web page at the 3XMM-DR5 catalog Website: http://xmm-catalog.irap.omp.eu/source/201125706010086.</footnote>

<footnote n="14">¹⁴ This model is appropriate for quasi-circular orbits; see http://www.atnf.csiro.au/people/pulsar/tempo/ref_man_seCTIONS/binary.txt.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Position of XB091D and its host globular cluster B091D in the Andromeda galaxy in several scales, using colored images from the Digitized Sky Survey and HST ACS. On the zoomed inset in the bottom-left corner, we overplot the 95% confidence level X-ray positional uncertainties from the *Chandra* (red circle) and *XMM-Newton* (blue circle) catalogs.</figcaption></figure>

We note that when correcting photon arrival times using the orbital solution of Esposito et al. (2016), we were only able to get sufficient pulsation significance in three observations (0112570101, 0650560301, 0690600401), whereas in all remaining data sets the pulsating signal was not detected because it was smeared. At the same time, the orbital solution obtained in this study allows pulsations to be significantly detected in all observations listed in Table 2. Although these two orbital parameter estimates agree on the order of magnitude with each other, the solution presented here is more precise.

This analysis of combined data sets reveals a spin-up trend observed in XB091D, which is otherwise hard to detect in individual observations; see Figure 3. Despite at least one probable $\dot{P}$ sign change (note, e.g., a period increase in the next to last data set comprising observations from 2012 June to August), it is likely that the neutron star spins up. On average, the spin-up rate amounts to $\dot{P} \approx -5.7 \times 10^{-13} \text{ s s}^{-1}$ if we consider all nine period estimates obtained, or $\dot{P} \approx -7.1 \times 10^{-13} \text{ s s}^{-1}$ if we reject the most recent data set, which has the largest period uncertainty.

## 5. Spectral Analysis

We considered all imaging observations in which XB091D was in the field of view (fov) of the EPIC instruments. Data reduction was performed with the standard *epproc* and *emproc* pipelines from the *XMM-Newton* Science Analysis Software (SAS v.14) and using the latest calibration files. From

the raw event lists, we selected high energy photons (>10 keV for the pn and between 10 and 12 keV for the MOS) and built light curves with 50 s binning, in order to identify and filter out time intervals affected by high particle background. Threshold count rates were set at 0.4 ct s⁻¹ and 0.35 ct s⁻¹ for the pn and MOS, respectively. We then used the interactive *xmmselect* task from the SAS to extract source and background spectra from the filtered event lists, starting with ObsID 0690600401, which has the largest counts statistics.

In this observation, the source lies in thefov of all three EPIC cameras. We extracted source counts from a circular region of 25" radius, while for the background we selected events from a nearby region free of sources within the same CCD. We applied the (`FLAG==0`) & (`PATTERN==0`) filtering options during the selection of pn events, in order to have the spectrum of the highest possible quality. For MOS counts, we used the standard filtering flags (`#XMMEA_EM`) and (`PATTERN<=12`). The extracted spectra were then re-binned in order to have at least 40 and 25 counts per energy bin for the pn and MOS, respectively.

Spectral fitting was performed with XSPEC v12.9 (Arnaud 1996). The spectrum of XB091D can be described by an absorbed power law with exponential cutoff, *wabs* (`cflux*cutoffpl`) in XSPEC. The spectrum is hard, with a photon index $\Gamma = 0.20 \pm 0.5$, a cutoff energy $E_{\text{cut}} = 4.6 \pm 0.4$ keV, and low absorption, $n_H = 3.79 \times 10^{20} \text{ cm}^{-2}$, obtained from PN spectrum alone. We kept this value fixed in the simultaneous pn+MOS fit. This best-fit column density
<table><caption>**Table 2**
XMM-Newton Observations Used for Timing Analysis in this Study and Its Main Results</caption><thead><tr><th>ObsID</th><th>Obs. start<br>UTC</th><th>Exposure<br>(s)</th><th>P<sub>spin</sub><br>(s)</th><th>S/dof</th><th>Pulse depth</th></tr></thead><tbody><tr><td>0112570101</td><td>2002 Jan 06 18:45:45</td><td>64317</td><td>1.203898(49)</td><td>8.7</td><td>0.33 ± 0.07</td></tr><tr><td>0405320701</td><td>2006 Dec 31 14:24:50</td><td>15918</td><td>1.203731(11)</td><td>4.1</td><td>0.28 ± 0.08</td></tr><tr><td>0405320901</td><td>2007 Feb 05 03:44:24</td><td>16914</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0505720201</td><td>2007 Dec 29 13:42:13</td><td>27541</td><td>1.203738(10)</td><td>7.3</td><td>0.20 ± 0.05</td></tr><tr><td>0505720301</td><td>2008 Jan 08 07:01:05</td><td>27219</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0505720401</td><td>2008 Jan 18 15:11:47</td><td>22817</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0505720501</td><td>2008 Jan 27 22:28:21</td><td>21818</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0551690201</td><td>2008 Dec 30 03:27:52</td><td>21916</td><td>1.203662(8)</td><td>5.0</td><td>0.23 ± 0.05</td></tr><tr><td>0551690301</td><td>2009 Jan 09 06:19:54</td><td>21918</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0551690501</td><td>2009 Jan 27 07:23:03</td><td>21912</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0551690601</td><td>2009 Feb 04 13:21:03</td><td>26917</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0600660201</td><td>2009 Dec 28 12:42:54</td><td>18820</td><td>1.203675(13)</td><td>4.2</td><td>0.25 ± 0.08</td></tr><tr><td>0600660501</td><td>2010 Jan 25 02:39:14</td><td>19715</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0650560301</td><td>2011 Jan 04 18:10:16</td><td>33415</td><td>1.203651(8)</td><td>9.1</td><td>0.31 ± 0.06</td></tr><tr><td>0650560401</td><td>2011 Jan 15 00:16:57</td><td>24316</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0650560601</td><td>2011 Feb 03 23:58:12</td><td>23918</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0674210201</td><td>2011 Dec 28 01:07:36</td><td>19034</td><td>1.203634(13)</td><td>4.9</td><td>0.35 ± 0.09</td></tr><tr><td>0674210301</td><td>2012 Jan 07 02:47:01</td><td>15433</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0674210401</td><td>2012 Jan 15 15:00:38</td><td>19916</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0674210501</td><td>2012 Jan 21 12:22:03</td><td>17317</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0690600401</td><td>2012 Jun 26 06:29:43</td><td>122355</td><td>1.203698(4)</td><td>20.6</td><td>0.23 ± 0.03</td></tr><tr><td>0700380501</td><td>2012 Jul 28 15:16:27</td><td>11914</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0700380601</td><td>2012 Aug 08 23:08:08</td><td>23916</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0701981201</td><td>2013 Feb 08 22:19:55</td><td>23918</td><td>1.20373(13)</td><td>7.0</td><td>0.45 ± 0.14</td></tr></tbody></table>

Note. Available observations were split into nine blocks which were analyzed with the assumption that the neutron star's spin period $P_{spin}$ does not change much within them. $S/dof$ is the statistical significance of the obtained solution as defined in Appendix A.1. Pulse depth (also referred to as pulse fraction) is defined as quantity A there as well. Uncertainties represent $1\sigma$ confidence interval.

is in broad agreement with the expected value from Dickey & Lockman (1990) in the direction of M31. The reduced chi-square of the fit is 1.16 for 310 degrees of freedom. We present a plot of the best-fit folded spectrum and model for ObsID 0690600401 in Figure 4.

We then analyzed all other observations, discarding those data sets where the source extraction region overlapped with CCD gaps or columns of bad pixels, and applying less stringent filtering options, (FLAG==0)&&(PATTERN<=4). In any case, given the lower count statistics, for all but ObsID 0700380601 a simpler absorbed power-law model, wabs(cflux*powerlaw), is sufficient to fit the spectra adequately. We used the multiplicative component cflux, which returns the source's flux directly as a parameter of the fit, to characterize variations with the epoch of the unabsorbed 0.3–10 keV band luminosity, and of the hardness ratio HR, where HR = Flux<sub>5–10 keV</sub>/Flux<sub>0.3–5 keV</sub>, assuming a distance to M31 of 752 kpc (Riess et al. 2012) and isotropic emission (Figure 5). Intriguingly, the shape of the spectrum seems to be related to the luminosity of the source, where the harder spectra occur at higher luminosities; see Figure 6.

Trudolyubov & Friedhorsky (2004) model the source X-ray spectrum with a hard power law. Shaw Greening et al. (2009) similarly find its photon index to be 0.9 ± 0.1 and interpret it as the high-mass X-ray binary (HMXB) nature of the source. Esposito et al. (2016) note that for some observations the absorbed power-law model does not provide satisfactory results and fit a blackbody + power-law model, as well as a cutoff

power law, which agree with our results within the uncertainties. They also note a similar “harder-brighter” correlation between the source luminosity and its spectral shape.

## 6. Discussion

We first estimate the chance association probability of this X-ray source with the GC B091D projection without being physically associated, as Esposito et al. (2016) favor an interpretation of the pulsar not associated with the GC. The association was first suggested by Supper et al. (2001) based on ROSAT data and confirmed by Trudolyubov & Friedhorsky (2004) when the first XMM-Newton and Chandra data became available. These works, however, do not provide an assessment of the probability of the false association.

XB091D is included in 3XMM-DR5 with R.A. = 00:43:01.478, decl. = +41:30:16.94 (J2000) and a $1\sigma$ positional uncertainty of 0.′13. These coordinates are computed as a weighted average of all individual detections, most of them being largely off-axis. The coordinates from ObsID 0690600401, the one with the smallest off-axis angle (1′.82), are R. A. = 00:43:01.440, decl. = +41:30:17.30 (J2000), or 0.′.56 away from the averaged coordinates, with a $1\sigma$ positional uncertainty of 0.′.19. Chandra Source Catalog Release 1.1 (Evans et al. 2010) includes this source as CXO J004301.4+413016 with coordinates R.A. = 00:43:01.469, decl. = +41:30:16.80 (J2000), and the major semi-axis of the 95% confidence level error ellipse of
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Recovered pulse profiles of XB091D obtained from the search for coherent pulsations using the regularized *S* statistics in nine combined data sets listed in Table 2. The shift along the *Y* axis is arbitrary and is made for clarity.</figcaption></figure>

1".19. The Bologna catalog gives the B091D center coordinates as R.A. = 00:43:01.446, decl. = +41:30:17.15 (J2000), which is consistent with the cluster center we find from the *HST* image within 0."15 of its astrometric calibration uncertainty. Distances between the cluster center and the X-ray source coordinates are 0."41 for the averaged *XMM-Newton* position, 0."17 for the ObsID 0690600401 position, and 0."44 for the *Chandra* position. In our association probability calculations, we therefore assume that the X-ray source is located within 1" from the cluster center.

There are 38 GCs (including candidates) between a 14' and 16' projected distance from M31's center in the Bologna catalog, which yields the density of 6.0 × 10⁻⁵ GCs per square arcsec at this angular distance from the galaxy center. This gives a low probability (1.2 × 10⁻²) of having a GC within 1" of one of the 70 brightest X-ray sources from the 3XMM-DR5

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Evolution of neutron star spin period. Spin period estimates obtained from individual observations are displayed in light color. They fluctuate strongly; however, a typical uncertainty of an individual measurement (plotted separately on the left) is large due to poor photon statistics. What might be interpreted as a frequent change of *P* is, in reality, driven by these uncertainties. For this reason, we use the combined data sets listed in Table 2 to estimate better local period solutions. These are indicated by black dots in the plot. Their uncertainties are computed as per Appendix A.1 and illustrate the better quality of the combined data sets analysis. Note that in some cases the spin period estimates are not available from individual observations; however, combining these data sets with adjacent ones often improves the pulse statistics. The dashed gray line is a linear fit to the period estimates from the combined data sets. The dashed cyan line is the same but discards the last period estimate due to its large uncertainty.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Best-fit X-ray spectrum and folded model of XB091D during ObsID 0690600401. pn data are black, MOS1 are green, and MOS2 are red. See the main text for details.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (Top) Unabsorbed 0.3–10 keV luminosity of XB091D, in units of $10^{38}$ erg s$^{-1}$, vs. the epoch of the observation. (Bottom) Ratio of the unabsorbed fluxes estimated in the 5–10 keV and 0.3–5 keV bands, respectively. Error bars correspond to $1\sigma$ confidence limits.</figcaption></figure>

catalog in M31. In fact, only 4 of the 70 brightest X-ray sources lie in this annulus. Therefore, their density of $6.3 \times 10^{-6}$ per square arcsec is translated to the probability of only $7.5 \times 10^{-4}$ of having a bright X-ray source within 1" of any of the GCs. It is therefore highly likely that this X-ray source belongs to the B091D GC. Given the probability of coincidence, we assume below that the source belongs to the GC, and do not discuss alternative interpretations.

The detection of pulsations at ~0.83 Hz secures the identification of the source with a spinning neutron star. This is the first object of this class identified in the Andromeda Galaxy, and one of the most distant pulsars observed to date. Also, given that XB091D resides inside the GC, it spins slower than any known sources in a GC (see the Introduction) and about 10 times slower than the slowest accreting pulsar in a GC known previously, the ~11 Hz IGR J17480-2466 in Terzan 5 (Papitto et al. 2011). Its orbital period of 30.5 hr is also the longest known for accreting GC binaries (see Table 5 in Bahramian et al. 2014).

Given the maximum observed X-ray luminosity of the system, $L_X = 1.2 \times 10^{38}$ erg s$^{-1}$ (see Figure 5), which is close to Eddington, one can determine the corresponding mass accretion rate assuming that bolometric luminosity does not greatly exceed the X-ray luminosity: $\dot{M} = \frac{L_X R_{NS}}{GM_{NS}} = 1.0 \times 10^{-8} M_{\odot} \text{yr}^{-1}$. We note that the observed X-ray variability of a factor ~3–4 (see Figure 5) excludes significant contribution from other unresolved X-ray sources residing in the same host GC. We cannot completely rule out the possibility that we observe two

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Hardness ratio vs. the unabsorbed 0.3–10 keV luminosity of XB091D. Error bars correspond to $1\sigma$ confidence limits.</figcaption></figure>

superimposed X-ray sources residing in the same GC, although the probability of this coincidence is relatively low: in our Galaxy there is only one GC known to host two bright persistent LMXBs, M15 (White & Angelini 2001). In case of the superposition, we can only refer to the pulsed part of the flux (typically 20%; see Table 2) as originating from this X-ray source. This would adjust our calculation of the accretion rate and the magnetic field below by a small factor, without affecting the important conclusions.

For a neutron star spinning with a 1.2 s period, the corotation radius is $R_C = (GMP^2/4\pi^2)^{1/3} \simeq 1890 M_{1.4}^{1/3}$ km, where $M_{1.4}$ is the neutron star mass in units of $1.4 M_{\odot}$.

We employ the same approach as in Papitto et al. (2011) to estimate the lower and upper limits of the magnetic field of the neutron star given that it persistently accretes matter from the companion. It is reasonable to assume in this case that the accretion disk is truncated at the radius $R_{\text{in}}$ between the radius of the neutron star $R_{\text{NS}}$ and the corotation radius $R_C$: $R_{\text{NS}} < R_{\text{in}} \lesssim R_C = 1890 M_{1.4}^{1/3}$. The inner disk radius can be approximately defined (e.g., Ghosh & Lamb 1978) from the balance between the kinetic energy of the falling gas and the magnetic energy of the neutron star magnetosphere, $R_{\text{in}} \simeq 160 M_{1.4}^{1/7} R_6^{-2/7} L_{37}^{-2/7} \mu_{28}^{4/7}$, where $R_6$ is the neutron star radius in units of 10 km, $L_{37}$ is the accretion luminosity in units of $10^{37}$ erg s$^{-1}$, and $\mu_{28}$ is the magnetic dipole moment of the neutron star in units of $10^{28}$ G cm$^3$. After trivial transformations, we obtain $0.008 L_{37\text{max}}^{1/2} < \mu_{28} \lesssim 75 L_{37\text{min}}^{1/2}$. We substitute the minimum ($L_{37\text{min}} = 3.5$) and maximum ($L_{37\text{max}} = 12.0$) observed X-ray luminosity into the resulting expression and obtain $2.7 \times 10^8 \text{ G} < B \lesssim 1.4 \times 10^{12} \text{ G}$. This upper limit corresponds to the case of equilibrium rotation of the neutron star when the accretion disk is truncated at the corotation radius and the falling matter does not transfer any angular momentum to the neutron star. If the system is indeed in equilibrium, the derivative of the neutron star spin period $\dot{P}$ would be
fluctuating in sign with small typical values. This is indeed what happens if one analyzes individual observations of XB091D (see Figure 3, light gray). Assuming coherent rotation between adjacent observations separated by a few months, we were able to significantly improve the statistics and constrain the period to a much better precision, which allowed us to detect a constant average spin-up of the neutron star.

We find an average spin-up rate $\dot{\nu} \approx 4.0-5.0 \times 10^{-13}$ Hz s$^{-1}$, only a factor of two to three less than that of IGR J17480-2446, which possesses a spin frequency derivative of $\dot{\nu} \approx 1.2 \times 10^{-12}$ Hz s$^{-1}$ (Cavecchi et al. 2011; Patruno et al. 2012). Such a significant average spin-up rate, if confirmed in future observations, may indicate that XB091D is not in equilibrium rotation, and in this case its magnetic field estimate should be lower than the equilibrium value of $B \simeq 1.4 \times 10^{12}$ G. IGR J17480-2446's average spin-up rate agrees well with the accretion of disk matter angular momentum given the observed luminosity $2-7 \times 10^{37}$ erg s$^{-1}$ and independent estimates of the inner radius of the accretion disk from the quasi-periodic oscillations (Papitto et al. 2011). From simple considerations, a torque exerted on the neutron star by accreted material moving in a Keplerian disk is $N = 2\pi I\dot{\nu} = \dot{M} \sqrt{GM R_{\text{in}}}$, where I is the neutron star moment of inertia usually assumed to be $10^{45}$ g cm$^2$. Therefore, we could roughly estimate the inner disk radius to be $R_{\text{in}} \simeq 30$ km and the corresponding magnetic field to be $B \simeq 5 \times 10^8$ G computed for the minimum observed accretion rate $\dot{M} = 2.6 \times 10^{-9} M_{\odot}$ yr$^{-1}$. This simple computation yields results that are very similar to the ones for IGR J17480-2446: $R_{\text{in}} \lesssim 20$ km, $B \simeq 7 \times 10^8$ G (Papitto et al. 2011). The observed average spin-up rate is not compatible with this simple model for the maximum observed accretion rate $\dot{M} = 1.0 \times 10^{-8} M_{\odot}$ yr$^{-1}$ as it gives $R_{\text{in}}$ below 10 km.

It is not clear to what extent this $\dot{P}$ estimate is influenced by the pulsar's acceleration along the line of sight $a_l$ in the GC gravitational field, which adds a term $a_l/c$ to $\dot{P}/P$. However, the largest known absolute value of $\dot{P}/P$ from pulsars in GCs (measured in B1718-19 in NGC 6342) amounts to $\dot{P}/P \approx 1 \times 10^{-15}$ s$^{-1}$ and can be used as the upper limit of $a_l/c$ for pulsars in GCs. We therefore do not expect it to influence our spin-up rate estimate by much.

On the other hand, many accreting X-ray pulsars (including LMXBs) are known to undergo lengthy episodes of spin ups and spin downs, at the same time staying close to their equilibrium rotation period (e.g., Bildsten et al. 1997; Chakrabarty et al. 1997; Perna et al. 2006). If confirmed, the penultimate data set in Figure 3 indicates at least one switch from the spin-up to the spin-down regime. We therefore cannot exclude the equilibrium rotation of the neutron star in XB091D and hence the larger value of its magnetic field $B \simeq 1.4 \times 10^{12}$ G.

XB091D exhibits rather small luminosity changes during the observed 11 year span (see Figure 5), which increases to 12 years if we consider an ObsID 0112570601 from the end of 2000 at a very large off-axis angle, which however shows a very similar flux in the 3XMM-DR5 catalog. The hard X-ray spectrum (see Figure 4) of the source is typical for LMXBs; however, it cannot be used alone to discard other possible interpretations, such as an HMXB nature. Keeping in mind the association with a GC, this evidence allows us to classify the source as a persistent LMXB. Assuming 100% duty cycle and persistent accretion at the observed rate, it becomes possible to estimate the time to spin up XB091D to a millisecond period, $P/\dot{P} \simeq 50-100 \times 10^3$ years. Extrapolating to the past, it

would take XB091D only 1 Myr or even less to spin up from a 10 s period to its observed value.

## 6.1. Properties of the Host Globular Cluster B091D and Comparison with Galactic Globular Clusters

In order to compare B091D to Terzan 5, the host cluster of the slowest known accreting X-ray pulsar in a GC (Papitto et al. 2011), and to several other massive GCs in the Milky Way hosting dense stellar cores, we used archival imaging and measurements of stellar velocity dispersion and stellar populations of GCs from the literature. Knowing the cluster structural parameters (e.g., the central stellar density $\rho_0$ and the core radius $r_c$) and central stellar velocity dispersions $\sigma_0$ in cluster cores, we can directly compare the numbers of LMXB formation events by tidal capture using the relation for the rate of encounters from Verbunt & Hut (1987):

$$ \Gamma \propto \rho_0^2 r_c^3 / \sigma_0. \quad (1) $$

Terzan 5 is a massive ($M = 2 \times 10^6 M_{\odot}$) very metal-rich cluster ([Fe/H] $\simeq 0$ dex) hosting a rich population of X-ray binaries. It possesses multiple stellar populations (Ferraro et al. 2009) so it is believed to be a nucleus of a dwarf galaxy heavily stripped by the Milky Way. It has a small core radius $r_c = 0.16$ pc, a very high central stellar density $\rho_0 = 1.58 \times 10^6 M_{\odot}$ pc$^{-3}$ (Prager et al. 2016, from the long-term timing observations of pulsars population), and a tidal radius $r_t \approx 6.7$ pc (Lanzoni et al. 2010). Its central velocity dispersion is estimated to be $\sigma_0 = 12.7$ km s$^{-1}$ (Gnedin et al. 2002). Terzan 5 has one of the densest known stellar cores among GCs in our Galaxy (Cohn et al. 2002).

B091D was included in the sample of M31 GCs with stellar velocity dispersions measured from high-resolution optical spectra (Strader et al. 2011). However, their dynamical models are based on ground-based images and therefore, might be inaccurate.

We downloaded high-resolution Hubble Space Telescope optical images for B091D from the Barbara A. Mikulski Archive for Space Telescopes<fnref n="15" /> obtained in the framework of the HST GO program 10273 “Accurately Mapping M31’s Microlensing Population” (P.I.: A. Crotts). The cluster is located close to the edge of thefov on two single exposures in the F555W (exposure time $t_{exp} = 151$ s) and F814W ($t_{exp} = 457$ s) filters obtained with the HST Advanced Camera for Surveys Wide Field Camera. Several pixels in the F814W image close to the cluster core are saturated. We generated HST point-spread functions in the two filters using the TINYTIM software (Krist et al. 2011) and then ran the GALFIT two-dimensional image-fitting code (Peng et al. 2002) and fitted King (1966) profiles into the B091D images, masking saturated pixels in the F814W image.

We obtained the following parameters for the F555W profile: core radius $r_c = 0.1110 = 0.42$ pc, truncation radius $r_t = 7.1 = 27$ pc, central surface brightness $\mu_{0.555} = 13.61$ mag arcsec$^{-2}$ or $I_0 = 1.31 \times 10^5 L_{\odot}$ pc$^{-2}$ which corresponds to $\rho_0 \approx 8 \times 10^5 M_{\odot}$ pc$^{-3}$, and ellipticity $e = 0.92$. Uncertainties on the structural parameters are on the order of 5%–7%. The parameters obtained from fitting the F814W image ($r_c = 0.115$, $r_t = 6.3$, $\mu_{0.814} = 13.44$ mag arcsec$^{-2}$, $e = 0.89$) are consistent within uncertainties with those obtained from F555W, except for the core radius, which is

<footnote n="15">¹⁵ http://mast.stsci.edu/</footnote>
probably affected by saturated pixels. We note that Agar & Barmby (2013) used the same data sets and performed a similar analysis. However, for some reason, their values of the central surface brightness and core radius derived from the F555W image differ from the F814W one by almost three orders of magnitude and look quite unrealistic. The latter set of parameters (F814W) agrees reasonably well with our estimates. We note that our $r_c$ value is somewhat smaller than the core radius $r_0 = 0.54$ pc reported by Peacock et al. (2010) and obtained from the analysis of ground-based near-infrared K-band images.

Following Richstone & Tremaine (1986) and converting into proper units, we estimate the V-band dynamical mass-to-light ratio of B091D to be $(M/L)_V \approx 333\sigma_0^2/(r_c I_0)$, where $\sigma_0$ is a central projected stellar velocity dispersion in km s$^{-1}$, $r_c$ is a core radius in pc, and $I_0$ is a central surface brightness in $L_\odot$ pc$^{-2}$. The aperture correction for the observed value $\sigma = 18.6 \pm 1.0$ km s$^{-1}$ (Strader et al. 2011) obtained by the integration of the King model yields $\sigma_0 = 21.0 \pm 1.3$ km s$^{-1}$, hence $(M/L)_{V,dyn} = 2.6 \pm 0.4(M_\odot/L_\odot)_V$ or $M_{dyn} = (9.6 \pm 1.5) \times 10^5 M_\odot$, assuming $V = 15.39$ mag.

Substituting these values into Equation (1), we estimate the ratio of the stellar encounter rates in Terzan 5 and B091D to be $\Gamma_{Terzan 5}/\Gamma_{B091D} \approx 0.34$, which means that LMXB formation by capture in the core of B091D is 2.9 times more likely than in Terzan 5, one of the densest GCs in our Galaxy, which also possesses the richest population of X-ray sources observed in a GC (Heinke et al. 2006). As the encounter number $\Gamma$ was shown to correlate with the number of X-ray binaries in GCs (Pooley et al. 2003), we can therefore expect B091D to be a quite prolific GC having an X-ray binary population similar to or larger than that of Terzan 5 and continuously forming new systems at the present epoch.

Caldwell et al. (2011) reported an old stellar population ($t \approx 12$ Gyr) and a metallicity [Fe/H] = -0.70 dex for B091D, making it representative of “red” metal-rich GCs. For the Kroupa (2002) stellar initial mass function, these parameters correspond to the stellar mass-to-light ratio $(M/L)_{*V} = 2.6(M_\odot/L_\odot)$. This remarkable agreement with the dynamical $(M/L)_{V,dyn}$ suggests that B091D did not experience stellar mass loss due to tidal stripping by M31 after dynamical relaxation and stellar mass segregation if it had been born with the Kroupa IMF.

Assuming a simple spherically symmetric model, the estimated two-body relaxation timescale of B091D is $t_{relax} \approx 1$ Gyr (Spitzer & Hart 1971; Spitzer 1987). Therefore, it is not expected to have undergone core collapse because it occurs at $t_{cc} \approx 17t_{relax}$ (Baumgardt et al. 2002). The dense stellar core in B091D probably formed over 10 Gyr ago and is not a result of its secular evolution.

We use the recent dynamical analysis for a large sample of Galactic GCs by Baumgardt (2017) and additional data from the Harris catalog (Harris 1996), 2010 edition,<fnref n="16" /> to compare B091D and similar clusters from the Milky Way (see Table 3). It turns out that none of the Galactic GCs reaches the central velocity dispersion $\sigma_0 = 21$ km s$^{-1}$ observed in B091D. The two closest analogs are NGC 6388 and NGC 6441. Interestingly, similarly to B091D, they also both harbor metal-rich stars (-0.55 and -0.46 dex), are massive (masses slightly over $10^6 M_\odot$), and possess dense stellar cores with central densities

<table><caption>**Table 3**
Structural Parameters of B091D and Similar Globular Clusters in the Milky Way, Their Stellar Encounter Rate $\Gamma$, and Their Stellar Encounter Rate for a Single Binary $\gamma$</caption><thead><tr><th>Name</th><th>log ρ<sub>0</sub><sup>a</sup><br>(L<sub>☉</sub>/pc<sup>3</sup>)</th><th>M/L<sup>b</sup><br>(☉)</th><th>r<sub>c</sub><sup>a</sup><br>(pc)</th><th>σ<sub>0</sub><sup>b</sup><br>(km s<sup>-1</sup>)</th><th>Γ</th><th>γ</th></tr></thead><tbody><tr><td>M4</td><td>3.64</td><td>1.70</td><td>0.72</td><td>4.5</td><td>1</td><td>1</td></tr><tr><td>Terzan 5</td><td>5.78<sup>c</sup></td><td>2.6<sup>c</sup></td><td>0.16<sup>c</sup></td><td>13<sup>d</sup></td><td>170</td><td>73</td></tr><tr><td>NGC6388</td><td>5.37</td><td>2.11</td><td>0.38</td><td>15</td><td>200</td><td>20</td></tr><tr><td>ω Cen</td><td>3.15</td><td>2.54</td><td>3.45</td><td>17</td><td>6.7</td><td>0.12</td></tr><tr><td>NGC6441</td><td>5.26</td><td>2.30</td><td>0.51</td><td>17</td><td>300</td><td>15</td></tr><tr><td>M54</td><td>4.69</td><td>2.18</td><td>0.62</td><td>19</td><td>30</td><td>3.4</td></tr><tr><td>B091D</td><td>5.49<sup>c</sup></td><td>2.6<sup>e</sup></td><td>0.42<sup>e</sup></td><td>21<sup>e</sup></td><td>490</td><td>23</td></tr></tbody></table>

**Notes.** $\omega$ Cen does not have a compact core but is included here because of its high central velocity dispersion. M4 is included only as the $\Gamma$ and $\gamma$ unit scale. $\Gamma$ is computed as $(\rho_0 M/L)^{2/3}/\sigma_0$, and $\gamma$ is computed as $(\rho_0 M/L)/\sigma_0$, where $\rho_0$ is a central volume luminosity density; the result is rounded to two significant figures. The list is sorted by central velocity dispersion.

a Harris (1996), 2010 edition, unless noted otherwise.

b Baumgardt (2017), unless noted otherwise.

c Prager et al. (2016), M/L in the center.

d Gnedin et al. (2002), King model estimate assuming $M/L = 3$.

e This study.

around $4-5 \times 10^5 M_\odot$ pc$^{-3}$, which were probably reached primordially because both clusters have expected core-collapse times exceeding the Hubble time. Two other clusters, although having much lower stellar metallicities, with central velocity dispersions approaching 20 km s$^{-1}$, analyzed by Baumgardt (2017) are $\omega$ Cen and M54. The latter one also possesses a compact dense stellar core similar to B091D. A common characteristic of these four objects and B091D is their high dynamical mass-to-light ratios, which suggest that their present day stellar mass functions are compatible with the Kroupa IMF for B091D, NGC 6388, and NGC 6441, and is even slightly steeper at low masses for M54 and $\omega$ Cen. Conversely, the vast majority of other Galactic GCs are substantially less massive. Their low-mass stars, therefore, have been depleted by tidal stripping because they are globally mass segregated, so that low-mass stars migrated to the outskirts, and hence their dynamical mass-to-light ratios are low compared to what one would expect for the Kroupa IMF.<fnref n="17" />

It is agreed in the community that the well-studied $\omega$ Cen (Gnedin et al. 2002; Bekki & Freeman 2003) and M54 (Siegel et al. 2007) are likely the tidally stripped nuclei of dwarf galaxies; M54 is thought to be residing at the core of the Sagittarius dwarf spheroidal galaxy. This makes them analogous to ultra-compact dwarf galaxies (UCDs; Drinkwater et al. 2003) observed in the nearby Virgo and Fornax clusters, which are known to have mass-to-light ratios consistent with the Kroupa IMF when the two-body relaxation timescales are long (Chilingarian et al. 2011). UCDs have been proven to be tidally stripped galactic nuclei by the recent discovery of central massive black holes (Seth et al. 2014) and by their stellar population properties (Francis et al. 2012), which are similar to much more massive, tidally stripped M32-like compact elliptical galaxies (Chilingarian et al. 2009;

<footnote n="17">¹⁷ It should be understood that the compact cores themselves have much shorter two-body relaxation times, typically tens of Myr, and they are probably mass segregated and dynamically evolved, as demonstrated by Heyl et al. (2015) based on a population of young white dwarfs in 47 Tuc.</footnote>

<footnote n="16">¹⁶ Electronic version available online at http://www.physics.mcmaster.ca/~harris/mwgc.dat.</footnote>
Chilingarian & Zolotukhin 2015). The only known stellar systems where central luminosity densities exceed $10^5 L_\odot pc^{-3}$ are nuclear star clusters in low- and intermediate-luminosity galaxies, cEs, some UCDs (Evstigneeva et al. 2008), and some GCs. Among ~160 Galactic GCs in the Harris catalog, only a handful possess $\rho_0 > 10^5 L_\odot pc^{-3}$ and are not marked as core collapsed. Three of them are Terzan 5, NGC 6388, and NGC 6441. In light of new structural and dynamical data (namely $\sigma_0, \rho_0, M/L$ ratios) presented in Baumgardt (2017) and also keeping in mind that stripped nuclei should contribute to the high end of the GC luminosity function (see e.g., Pfeffer et al. 2014), we suspect that NGC 6388 and NGC 6441 are, in fact, heavily tidally stripped galactic nuclei. Interestingly, these two clusters have extremely extended and peculiar horizontal branches, which can be explained only by very high helium abundances (see Bellini et al. 2013; Brown et al. 2016) that may be evidence of self-enrichment during the extended period of star formation common in galactic nuclei but uncommon in GCs; this makes them unique compared to other GCs (Brown et al. 2016). They are similar to B091D in terms of both stellar population and internal structure and dynamics. They, however, must have undergone a somewhat higher degree of tidal stripping than M54.

Keeping in mind that Terzan 5 is hosting another slow X-ray pulsar that is also suspected to be a stripped galaxy nucleus, this raises the question whether some specific conditions in nuclear star clusters favor the formation of LMXBs. It is possible that the very high stellar densities reached in galactic nuclei, similar to those in core-collapsed GCs but on an order of magnitude larger spatial scale (see Table 2 in Evstigneeva et al. 2008 for central surface brightness and size values of UCD cores), provide an effective formation channel for relativistic binaries via close encounters (see also Dabringhausen et al. 2012 on the overabundance of LMXBs in UCDs). We illustrate this by comparing the $\Gamma$ values of the discussed clusters in Table 3. NGC 6388, NGC 6441, and B091D have large central densities $\rho_0$ and at the same time large cores $r_c$, and therefore very high encounter rates $\Gamma$ that even exceed those of the densest cluster Terzan 5.

We also compare another dynamical parameter for GCs, the encounter rate for a single binary, $\gamma \propto \rho_0/\sigma_0$ (Verbunt 2003; Verbunt & Freire 2014). It is expected that a higher $\gamma$ indicates a higher rate of exchange encounters in a GC and therefore a higher observed frequency of exchange encounter products— such as isolated pulsars, slow young neutron stars, and other kinds of exotic objects believed to be formed by the disruption of X-ray binaries. The lifetime of a binary until the next encounter, which increases their chances to get disrupted or to exchange companion stars, is proportional to $1/\gamma$. We computed $\gamma$ values for B091D, Terzan 5, and the four clusters discussed above, which we think can all be classified as UCDs using the updated global mass-to-light ratios and central velocity dispersions from Baumgardt (2017) and structural parameters from the Harris catalog. We note that Baumgardt (2017) demonstrated (see their Figure 3) that the $M/L$ ratio in the cluster core is close to the global value while the radial $M/L$ profiles often exhibit a dip explained by mass segregation. In the units of the reference GC M4 from Verbunt & Freire (2014), $\gamma_{B091D} = 23 \gamma_{M4}$. This makes B091D similar to the top five GCs of our Galaxy using this parameter.

## 6.2. System Age and Formation Scenarios

Caldwell et al. (2011) estimated the age of the B091D GC to be 12 Gyr. Generally, there exist few evolutionary sequences of binaries formed more than 10 Gyr ago that start an accretion episode of required intensity at the present time and hence could explain the origin of XB091D. For instance, in Podsiadlowski et al. (2002), binaries with the initial mass of the secondary star $M_2$ between 1.0 and 1.2 $M_\odot$ and initial orbital period $P$ between 0.5 and 100 days exhibit accretion episodes of order $\simeq 100$ Myr in duration after $\simeq 10–12$ Gyr of evolution, reaching the peak accretion rate $\dot{M}_{\text{peak}} \simeq \text{few} \times 10^{-8} M_\odot \text{ yr}^{-1}$, with an average accretion rate $\langle M \rangle \simeq \text{few} \times 10^{-9} M_\odot \text{ yr}^{-1}$. This shows that scenario of the primordial origin of XB091D, when the system formed at the early epochs of its host GC around 12 Gyr ago and very recently started the accretion episode that we observe today, is not forbidden by the evolution theories of isolated binary systems.

A scenario where XB091D hosts a primordial neutron star must explain how a neutron star that formed 12 Gyr ago retained the most probable value of its current magnetic field $B \simeq 1.4 \times 10^{12}$ G, provided that neutron stars are thought to be born with magnetic fields $B \simeq 10^{13–10^{14}}$ G (Faucher-Giguère & Kaspi 2006). It is not clear whether or not non-accreting neutron stars preserve their magnetic field for a long time. Some models estimate it to decay with a characteristic time that spans from $\sim 10^6$ (e.g., Narayan & Ostriker 1990) to $\sim 10^8$ years (e.g., Bhattacharya et al. 1992). Given the lack of consensus about the magnetic field decay in neutron stars, we cannot rule out the primordial neutron star origin in XB091D. In this context, we note that Ivanova et al. (2008) find in their simulations that less than 10% of the primordial binaries in a given GC that survived a core-collapse supernova remain in the original system after 11 Gyr (see their Table 4). In this case, if XB091D hosts a primordial neutron star, it must have acquired the current secondary at a later epoch. Then, a probable channel to get a companion is a binary exchange (Verbunt & Freire 2014). This is supported in particular by the fact that exchange encounters favor wide binary systems like XB091D with orbital periods of more than 1 day (e.g., Verbunt 2003).

An alternative formation scenario that is capable of producing a neutron star in an old GC is an accretion-induced collapse (AIC), wherein a massive ($\simeq 1.2 M_\odot$) ONeMg white dwarf (WD) accretes matter from a companion until it reaches the Chandrasekhar limit $M = 1.44 M_\odot$. AIC is anticipated to be responsible for the population of slow isolated pulsars with high magnetic fields in GCs (e.g., Lyne et al. 1996; Breton et al. 2007) and to be the origin of some slow accreting X-ray pulsars in the field such as 4U 1626–67 (e.g., Yungelson et al. 2002), although it has never been observed directly. Ivanova et al. (2008) claim that, in a typical GC during the 9.5–12.5 Gyr production of LMXBs from AIC, the AIC is two to three times more efficient than any other dynamical formation channel, such as physical collisions, tidal captures, and binary exchanges. A cluster like Terzan 5 (and therefore very similar to B091D) is thought to produce via AIC $9.35 \pm 1.20$ LMXBs per Gyr at ages $11 \pm 1.5$ Gyr (Ivanova et al. 2008), which is not negligible even considering the short lifetime of such binaries. One characteristic property of a neutron star formed in an AIC event is its low mass, $M_{\text{NS}} \simeq 1.26 M_\odot$. We doubt however that it is possible to constrain the neutron star mass in this system to support or discard an AIC hypothesis. No optical identification or
spectroscopic observations of the companion star of a source in the GC at the distance of M31 is possible with the current generation of astronomical instrumentation. In the case of an AIC origin of the neutron star in XB091D, there are several reasons to suspect that we are currently observing an accretion episode powered by a new donor after a dynamical exchange event with the original binary that hosted AIC took place. During AIC, a white dwarf loses roughly 0.2 $M_\odot$ in the form of binding energy and probably some mass in a supernova shell, which makes the binary orbit wider, therefore detaching the binary and halting mass transfer. The time between AIC and resumption of mass transfer in the ultra-compact system with a white dwarf donor is $\simeq 10^8$ years (Verbunt et al. 1990), though it obviously strongly depends on the donor properties—for example, on the presence of magnetic braking which brings the secondary into contact with its Roche lobe, or the rate of the secondary star radius increase due to its nuclear evolution. On the other hand, the donor loses a significant fraction of its mass to a power AIC event so that it can be incapable of powering an intensive accretion episode, which we are currently observing. Prolonged epochs without accretion increase the chances for a binary exchange inside a GC with a high specific encounter rate $\gamma$.

The mass loss of the original donor to trigger the AIC event should be within the 0.2–0.3 $M_\odot$ range. From our orbital solution, we estimate the donor mass function to be 0.0160, in agreement with Esposito et al. (2016). This is translated to the minimum donor mass of $M_2 = 0.36 M_\odot$ in the case of an edge-on system with inclination angle $i = 90$ deg. The lack of X-ray eclipses in the longest observation, which covers all orbital periods, means we can constrain the donor mass a little further as the system’s inclination is then less than $\simeq 70$ deg: $M_2 > 0.38 M_\odot$. For a random distribution of inclination angles, one has a 90% a priori probability of observing a binary system at an angle $i > 26$ deg. For the observed mass function, this inclination corresponds to $M_2 = 1.04 M_\odot$. Therefore, the 90% confidence interval for the donor mass is $0.38 \le M_2 \le 1.04 M_\odot$. In fact, for a 12 Gyr old GC B091D, the main-sequence turnoff mass is $0.8 M_\odot$. All stars within B091D with mass $\gtrsim 1.0 M_\odot$ must have evolved to red giants and even become white dwarfs. Thus, the reasonable upper limit on the donor mass is $M_2 \simeq 0.8–0.9 M_\odot$, which corresponds to the low-mass subgiants and stars leaving the main sequence.

Using the Eggleton (1983) formula and considering limits for the mass ratio ($M_1 = 1.25 M_\odot$, $M_2 = 0.9 M_\odot$, and $M_1 = 2.0 M_\odot$, $M_2 = 0.4 M_\odot$), it is easy to estimate the size of the Roche lobe for the companion star, $1.64 \le R_{L2} \le 2.24 R_\odot$. Therefore, to fill its Roche lobe, the donor must be an evolved star that has recently left the main sequence and has had its radius increased to about $2 R_\odot$. The current high accretion rate can be understood as driven by the nuclear evolution of the binary with a rapidly increasing radius.

It is very unlikely that a low-mass star with $M_2 \simeq 0.4–0.5 M_\odot$ (e.g., the donor that powered the AIC event) can reach a $1.6 R_\odot$ radius. This could, however, be influenced by the X-ray irradiation of the donor surface, but its effect is not well-understood presently. For the most likely value of the donor mass $M_2 = 0.8 M_\odot$, the inclination of the orbital plane is $\simeq 30$ deg and the orbital separation is $6.27–6.96 R_\odot$ (0.029–0.032 au). This makes XB091D the widest known accreting binary system in a GC.

The pulsar recycling theory (Bhattacharya & van den Heuvel 1991) assumes that after the rotation-powered phase of a classical pulsar finishing with $P_{\text{spin}} \approx 1–10$ s, neutron stars can be spun up to $P_{\text{spin}} \approx 1–10$ ms by an accreting donor in binary system. XB091D fits very well into this picture being observed at the earliest stages of its accretion spin-up phase. Whereas a very similar pulsar, IGR J17480–2466, from Terzan 5 represents a mildly recycled system, XB091D is a missing example of a non-recycled neutron star, which nonetheless is accreting very intensively. The endpoint of the evolution of XB091D in a few million years is likely to be a millisecond radio pulsar in a wide ($\simeq$several days) orbit with a white dwarf companion. Such MSP systems are numerous, but their progenitors like IGR J17480–2466 and XB091D have relatively short lifetimes at high mass transfer rates, so very few of them are observed.

## 7. Conclusions

We report an independent detection of a luminous ($L_X = 3–12 \times 10^{37}$ erg s$^{-1}$) accreting X-ray pulsar in the Andromeda galaxy in the public data of 38 observations obtained by the XMM-Newton observatory between 2000 and 2013. In 13 observations, we detected 15% to 30% pulsed emission with a period of 1.2 s. Our analysis is fully reproducible online using the XMM-Newton photon database available at http://xmm-catalog.irap.omp.eu.

We demonstrate that this X-ray binary is associated with a massive 12 Gyr old GC B091D hosting a very dense stellar core and possessing a high stellar encounter rate. Therefore, this system is a very unusual example of a non-recycled pulsar intensively accreting. At 1.2 s, its neutron star spins 10 times slower than the former slowest known X-ray pulsar in GCs—a mildly recycled system, IGR J17480–2446, in Terzan 5.

From the X-ray timing analysis, we estimate the binary system orbital parameters, including its orbital period of 30.5 hr. By combining several adjacent data sets in order to increase the photon statistics, we obtain a phase-connected solution in nine extended periods of time over the 11 year baseline and detect an average neutron star’s spin-up rate of $\dot{P} \simeq -7.1 \times 10^{-13}$ s s$^{-1}$, which however has at least one episode of spin down. If the system is not in an equilibrium rotation and the spin up persists, from this number we estimate that the accretion onset happened less than 1 Myr ago, because it will only take $\approx 10^5$ years for this system to become a conventional millisecond pulsar.

The detected average spin-up rate matches that expected from the angular momentum accretion of a Keplerian disk with the accretion rate $\dot{M} \simeq 3 \times 10^{-9} M_\odot$ yr$^{-1}$ if the inner boundary of the disk is at $R_{\text{in}} \simeq 30$ km. If we assume that the neutron star in XB091D indeed is not in equilibrium rotation, we can estimate its magnetic field to be $B \simeq 5 \times 10^8$ G. The observed change to spin down, however, favors an equilibrium configuration with a larger value $B \simeq 1.4 \times 10^{12}$ G, which is also supported by the hard observed X-ray spectrum. From the orbital separation and the donor Roche lobe size, also keeping in mind that the system has been persistently accreting over the past 12 years, we conclude that the secondary must be a slightly evolved low-mass star with the mass close to the main-sequence turnoff for a 12 Gyr old GC $M_2 \simeq 0.8 M_\odot$.

Based on these properties, we cannot distinguish between a primordial neutron star origin and its formation at a later epoch, e.g., in the AIC event. However, in both cases, it is highly
unlikely that we observe an original binary where the neutron star was formed. In the AIC case, the system has likely experienced an exchange interaction and the neutron star captured a non-exhausted low-mass donor star. In the case of a primordial neutron star, it must have acquired a secondary star after it was formed. In both cases, the most likely process of getting a donor star is an exchange encounter and the donor started to overflow its Roche lobe very recently, less than 1 Myr ago. These formation scenarios are in line with the measured properties of XB091D and correspond to the expectations that follow from the global properties of the host GC, namely its high encounter rate for a single binary $\gamma$, a predicted indicator of the frequency of binary systems that form via exchange encounters. XB091D is the first accreting non-recycled X-ray pulsar, completing the picture of pulsar recycling.

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France. Part of the plots were produced using Veusz by Jeremy Sanders. The authors are grateful to citizen scientists M. Chernyshov, A. Sergeev, and A. Timirgazin for their help with the development of the XMM-Newton catalog Website, http://xmm-catalog.irap.omp.eu, used throughout this study. The authors thank N. Ivanova for the useful comments on the paper and H. Baumgardt for providing structural parameters for globular clusters and discussion regarding the globular cluster dynamics. I.Z. acknowledges the support by the Russian Scientific Foundation grant 14-50-00043 for the data processing and grant 14-12-00146 for the timing analysis. I.C. and I.Z. acknowledge the joint RFBR/CNRS grant 15-52-15050 supporting the Russian-French collaboration on the archival and Virtual Observatory research, the RFBR grant 15-32-21062, and the president of the Russian Federation grant MD-7355.2015.2 supporting the studies of globular clusters and compact stellar systems. The work of N.S. was supported by the French Space Agency CNES through the CNRS. M.B. was supported by the Sardinian Region through a fundamental research grant under Regional Law 7th. The authors also acknowledge support by Paris Astronomical Data Center (PADC). Part of the detection chain was adapted from the software library for X-ray timing, MaLTPyNT (Bachetti 2015).

# Appendix

## Improved Period Search and Error Analysis for the Epoch-folding Technique

### A.1. Error Analysis for the Epoch Folding of a Single Observation

Here we follow Leahy et al. (1983) and Leahy (1987) in order to perform the period search using the epoch-folding technique and make an additional step and estimate the period determination uncertainty from statistical considerations.

Let us consider a harmonic signal with a period $P$ over some constant background so that the pulse shape is expressed as

$$ 
\begin{aligned}
f(t) &= a + b \sin t \\
a &= N_{\gamma}/T,
\end{aligned} 
\qquad (2) 
$$

where $b$ is the pulse amplitude and $a$ is the background value that we estimate from the total number of photons $N_\gamma$ registered during the total exposure time $T$. Here we assume that the exposure-filling factor is 100%, i.e., no gaps took place during the exposure time due to, e.g., soft proton flares.

If we now perform the epoch folding with a slightly different period $P + \Delta P$, $\Delta P \ll P$, it will cause a phase shift of the last pulse:

$$ \Delta\phi = \frac{2\pi T \Delta P}{P^2}. \qquad (3) $$

Hereafter, we take the continuum limit and replace all sums with integrals. The discretization does not change the final results much because the reduction in the method sensitivity is only 3.3% for $n = 10$ bins per phase and 0.8% for the $n = 20$ bins we used here (Leahy 1987). Then, we can estimate the phase-smeared folded pulse shape as a function of $\Delta\phi$:

$$ 
\begin{aligned}
f_1(t, \Delta\phi) &= \frac{1}{\Delta\phi} \int_0^{\Delta\phi} (a + b \sin(t - \tau)) d\tau \\
&= a + \frac{2b}{\Delta\phi} \sin \frac{\Delta\phi}{2} \sin \left( t - \frac{\Delta\phi}{2} \right).
\end{aligned} 
\qquad (4) 
$$

Now we can compute the $S$ statistics (Leahy et al. 1983) as

$$ 
\begin{aligned}
S(\Delta\phi) &= \int_{0}^{2\pi} \frac{(f_1(t, \Delta\phi) - a)^2}{f_1(t, \Delta\phi)} dt \\
&\approx \frac{1}{a} \int_{0}^{2\pi} (f_1(t, \Delta\phi) - a)^2 dt \\
&= \frac{2\pi b^2}{a \Delta\phi^2} (1 - \cos \Delta\phi) = \frac{4\pi b^2}{a \Delta\phi^2} \sin^2 \frac{\Delta\phi}{2}.
\end{aligned} 
\qquad (5) 
$$

For the simplicity of the computation, here we assume that the pulse is shallow (e.g., $b \ll a$) and hence we take $1/a$ outside the integral. More general analytical calculation is bulky and does not change the result very much because it depends weakly on the $b/a$ ratio as $\sim \sqrt{1 - (b/a)^2}$.

Now keeping in mind that $S$ is in fact the $\chi^2$ statistics, we can estimate the period uncertainty by solving the equation $S(\Delta\phi) = S(0) - 1$. The Taylor expansion of Equation (5) to the fourth power on $\Delta\phi$ yields

$$ \frac{\pi b^2}{a} \left( 1 - \frac{\Delta\phi^2}{12} \right) = \frac{\pi b^2}{a} - 1. \qquad (6) $$

Solving it for $\Delta\phi$ and introducing the pulse depth $A = b/a$ (as in Leahy 1987),

$$ \Delta P = \sqrt{\frac{3}{\pi^3} \frac{P^2}{A \sqrt{N_\gamma T}}}. \qquad (7) $$

### A.2. Analytic Formulation for the Epoch Folding of Coherent Pulsations in Two Observations

Now let us consider a harmonic signal with a period $P$ observed in two observations with exposure times $T$ and $T/n$ (without the loss of generality, we take a real number $n \ge 1$). The second observation starts at the moment $mT$ ($m$ is a real number, $m \ge 1$). Then, the folded phase shape resulting from
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Examples of the S statistics for the period search (without the stochastic component) for one observation (black line) with a total exposure T shown as the black solid line and for two observations (red line) with exposure times T₁ = T and T₂ = T/4. with the second observations starting at the moment T₃ = 12T. Shaded areas show examples of confidence regions: for the red curve it contains two secondary peaks on either side of the primary one.</figcaption></figure>

the sum of the two observations becomes

$$f_1(t, \Delta\phi, m, n) = \frac{1}{\Delta\phi} \left( \int_0^{\Delta\phi} (a + b \sin(t - \tau)) d\tau + \int_{m\Delta\phi}^{(m+1/n)\Delta\phi} (a + b \sin(t - \tau)) d\tau \right). \quad (8)$$

Here the second integral has a multiplier (1/Δφ) rather than (n/Δφ) because it will contribute as (1/n) to the total pulse. Omitting bulky computations and trigonometric transformations, the S statistics computed in the shallow pulse approximation (b ≪ a; see above) becomes

$$S(\Delta\phi, m, n) = \frac{4\pi b^2}{a\Delta\phi^2} \left( \sin^2\frac{\Delta\phi}{2} + \sin^2\frac{\Delta\phi}{2n} + 2\sin\frac{\Delta\phi}{2}\sin\frac{\Delta\phi}{2n}\cos\left(\frac{\Delta\phi}{2} - \frac{\Delta\phi}{2n} - m\Delta\phi\right) \right). \quad (9)$$

This expression is non-negative for any Δφ, and it has several properties of interest for our analysis. It is virtually identical to the light pattern formed by the double-slit diffraction of a coherent source on slits of unequal widths. Adding the second data set separated from the first one in time introduces the modulation of the S statistics from a single observation (Equation (5)) with the amplitude (1 + 1/n)² and the high frequency m/2π (see Figure 7). Then, depending on the statistics defined by the total number of registered photons, the 1σ confidence region for the period P may either shrink into a single modulated peak so that the period uncertainty ΔP will

be given by an equation similar to Equation (7):

$$\Delta P_{\text{peak}} = \sqrt{\frac{3}{\pi^3} \frac{P^2 \sqrt{1 + 1/n}}{A \sqrt{N_{\gamma} T (m + 1/n)}}} \quad (10)$$

or cover several secondary peaks, in which case a single value of ΔP becomes meaningless because the real P value may reside in one of several secondary peaks in the confidence region. The multiple solution situation occurs when the phase shift Δφ that corresponds to the period uncertainty ΔP from Equation (5) (with T corrected to the total exposure time T<sub>tot</sub> = T(1 + 1/n)) exceeds the phase distance between secondary peaks, Δφ > 2π/m assuming m ≪ 1. Hence, we can estimate the number of secondary peaks n<sub>sec</sub> within the upper envelope of the S distribution by equating the phase shift to 2πn<sub>sec</sub>/m:

$$n_{\text{sec}} = 2 \left[ \sqrt{\frac{3}{\pi^3} \frac{m}{A\sqrt{a}}} \right]. \quad (11)$$

This quantity depends on the count rate a = N<sub>γ</sub>/T<sub>tot</sub> and not on the actual counts. This conclusion looks counterintuitive at first. However, it is trivially explained by the fact that we work in phase space. Therefore, when the pulse statistics grows because of the increased exposure time, the main peak becomes narrower in terms of ΔP, and if we keep the separation between observations constant, m will decrease. In the case of XB091D, a ≈ 0.1, and A ≈ 0.3, we will get ~3m secondary peaks on either side of the primary one if we combine two observations. Another important property of Equation (9) is that one can use the measured modulation amplitude of the S statistics for the period search in the two observations in order to estimate the possible period variation. If the observed modulation stays close to the value predicted by Equation (9)
<figure><img src="image_9.png" /><figcaption>Figure 8. Pulse profile recovered from the observation collected on 2002 January 06 (ObsID: 0112570101). (Left) The yellow and red lines are pulse shapes obtained with and without regularization, respectively. (Right) The green and blue lines represent the S statistics vs. spin period for the search with and without regularization, respectively.</figcaption></figure>

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

plus the stochastic component $S_{noise}$, which is a $\chi^2_{n-1}$ statistics with mean $(n-1)$ (Leahy 1987), we can conclude that the period change between two observations is undetectable.

All the computations presented above can be trivially generalized to the case of $N_{obs} > 2$. However, all expressions become very lengthy and difficult to understand. The envelope size for the confidence region that might still contain multiple narrow peaks ($1\sigma$) can be estimated as

$$ \Delta P_{\text{mult}} = \sqrt{\frac{3}{\pi^3}} \frac{P^2}{A\sqrt{N_{\gamma \text{tot}} T_{\text{tot}}}}, \quad (12) $$

where $T_{tot}$ indicates the total good time interval (GTI) and $N_{\gamma \text{tot}}$ is the total number of detected photons.

### A.3. Regularization of the Pulse Shape

We can improve the period search procedure by assuming the pulse shape to be smooth. This idea is similar to regularization techniques used in the image/signal reconstruction, and therefore we can exploit similar mathematical methods. As long as the period is determined by searching the maximum of the S statistics, we introduce the penalization factor that depends on the pulse shape $p(\phi)$ as

$$ S_{\text{reg}} = S \left( 1 + k \int_{0}^{2\pi} \left( \frac{d^2 p(\phi)}{d\phi^2} \right)^2 d\phi \right)^{1/2}. \quad (13) $$

Here k is a positive coefficient that defines the degree of regularization. Numerically, the second derivative is computed for a discrete pulse shape as $\frac{3}{2}(\text{box}(p, 3) - p)$, where box($p, 3$) denotes the boxcar smoothing of the pulse $p$ with a window of 3 pixels.

Applying the regularization to real data results in much smoother pulse shapes, while the periods always stay within the uncertainties predicted from the S statistics or the analytic formulations provided above (Equation (12)). We stress that we do not smooth the actual pulse shape but rather bias the period search statistically toward smoother shapes. In Figure 8 (left panel), we show an example of the recovered pulse shape from

the data set obtained in 2002 January (ObsID: 0112570101) with and without regularization using $k = (2/3)^2$. The right panel shows the S statistics without (blue) and with (green) regularization obtained from the folded epoch period search in 20 phase bins. Therefore, the S statistics includes a stochastic contribution with the variance $\langle S_{\text{rand}} \rangle = 19$. It is clear from the plot that the regularized solution lies within the $1\sigma$ confidence area, i.e., the period values are consistent within statistical uncertainties; however, the recovered pulse shape looks substantially smoother in the regularized case.

## References

* Agar, J. R. R., & Barmby, P. 2013, AJ, 146, 135

* Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Natur, 300, 728

* Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Sci, 324, 1411

* Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

* Bachetti, M. 2015, MaLTPyNT, Astrophysics Source Code Library, ascl:1502.021

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

* Bahramian, A., Heinke, C. O., Sivakoff, G. R., et al. 2014, ApJ, 780, 127

* Baumgardt, H. 2017, MNRAS, 464, 2174

* Baumgardt, H., Hut, P., & Heggie, D. C. 2002, MNRAS, 336, 1069

* Bekki, K., & Freeman, K. C. 2003, MNRAS, 346, L11

* Bellini, A., Piotto, G., Milone, A. P., et al. 2013, ApJ, 765, 32

* Bhattacharya, D., & van den Heuvel, E. P. J. 1991, PhR, 203, 1

* Bhattacharya, D., Wijers, R. A. M. J., Hartman, J. W., & Verbunt, F. 1992, A&A, 254, 198

* Bildsten, L., Chakrabarty, D., Chiu, J., et al. 1997, ApJS, 113, 367

* Breton, R. P., Roberts, M. S. E., Ransom, S. M., et al. 2007, ApJ, 661, 1073

* Brown, T. M., Cassisi, S., D'Antona, F., et al. 2016, ApJ, 822, 44

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

* Caldwell, N., Schiavon, R., Morrison, H., Rose, J. A., & Harding, P. 2011, AJ, 141, 61

* Cavecci, Y., Patruno, A., Haskell, B., et al. 2011, ApJL, 740, L8

* Chakrabarty, D., Bildsten, L., Grunsfeld, J. M., et al. 1997, ApJ, 474, 414

* Chilingarian, I., Cayatte, V., Revaz, Y., et al. 2009, Sci, 326, 1379

* Chilingarian, I., & Zolotukhin, I. 2015, Sci, 348, 418

* Chilingarian, I. V., Mieske, S., Hilker, M., & Infante, L. 2011, MNRAS, 412, 1627

* Cohn, H. N., Lugger, P. M., Grindlay, J. E., & Edmonds, P. D. 2002, ApJ, 571, 818

* Dabringhausen, J., Kroupa, P., Pflamm-Altenburg, J., & Mieske, S. 2012, ApJ, 747, 72
Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215

Drinkwater, M. J., Gregg, M. D., Hilker, M., et al. 2003, Natur, 423, 519

Eggleton, P. 1983, ApJ, 268, 368

Esposito, P., Israel, G. L., Belfiore, A., et al. 2016, MNRAS, 457, L5

Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS, 189, 37

Evrstigneeva, E. A., Drinkwater, M. J., Peng, C. Y., et al. 2008, AJ, 136, 461

Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15

Faucher-Giguère, C.-A., & Kaspi, V. M. 2006, ApJ, 643, 332

Ferraro, F. R., Dalessandro, E., Mucciarelli, A., et al. 2009, Natur, 462, 483

Francis, K. J., Drinkwater, M. J., Chilingarian, I. V., Bolt, A. M., & Firth, P. 2012, MNRAS, 425, 325

Galleti, S., Federici, L., Bellazzini, M., Fusi Pecci, F., & Macrina, S. 2004, A&A, 416, 917

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

Gnedin, O. Y., Zhao, H., Pringle, J. E., et al. 2002, ApJL, 568, L23

Harris, W. E. 1996, AJ, 112, 1487

Heinke, C. O., Wijnands, R., Cohn, H. N., et al. 2006, ApJ, 651, 1098

Hénon, M. 1961, AnAp, 24, 369

Heyl, J., Richer, H. B., Antolini, E., et al. 2015, ApJ, 804, 53

Hut, P., McMillan, S., Goodman, J., et al. 1992, PASP, 104, 981

Ivanova, N., Heinke, C. O., Rasio, F. A., Belczynski, K., & Fregeau, J. M. 2008, MNRAS, 386, 553

King, A., & Lasota, J.-P. 2016, arXiv:1601.03738

King, I. R. 1966, AJ, 71, 64

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Kroupa, P. 2002, Sci, 295, 82

Lanzoni, B., Ferraro, F. R., Dalessandro, E., et al. 2010, ApJ, 717, 653

Leahy, D. A. 1987, A&A, 180, 275

Leahy, D. A., Darbro, W., Elsner, R. F., et al. 1983, ApJ, 266, 160

Lyne, A. G., Biggs, J. D., Harrison, P. A., & Bailes, M. 1993, Natur, 361, 47

Lyne, A. G., Manchester, R. N., & D'Amico, N. 1996, ApJL, 460, L41

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

Michel, F. C. 1987, Natur, 329, 310

Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015, MNRAS, 447, 1847

Narayan, R., & Ostriker, J. P. 1990, ApJ, 352, 222

Olausen, S. A., & Kaspi, V. M. 2014, ApJS, 212, 6

Osborne, J. P., Borozdin, K. N., Trudolyubov, S. P., et al. 2001, A&A, 378, 800

Papitto, A., D'Ai, A., Motta, S., et al. 2011, A&A, 526, L3

Patruno, A., Alpar, M. A., van der Klis, M., & van den Heuvel, E. P. J. 2012, ApJ, 752, 33

Patruno, A., & Watts, A. L. 2012, arXiv:1206.2727

Peacock, M. B., Maccarone, T. J., Knigge, C., et al. 2010, MNRAS, 402, 803

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266

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

Pfeffer, J., Griffen, B. F., Baumgardt, H., & Hilker, M. 2014, MNRAS, 444, 3670

Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107

Pooley, D., Lewin, W. H. G., Anderson, S. F., et al. 2003, ApJL, 591, L131

Prager, B., Ransom, S., Freire, P., et al. 2016, arXiv:1612.04395

Radhakrishnan, V., & Srinivasan, G. 1982, CSci, 51, 1096

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

Richstone, D. O., & Tremaine, S. 1986, AJ, 92, 72

Riess, A. G., Fliri, J., & Valls-Gabaud, D. 2012, ApJ, 745, 156

Rosen, S. R., Webb, N. A., Watson, M. G., et al. 2015, arXiv:1504.07051

Seth, A. C., van den Bosch, R., Mieske, S., et al. 2014, Natur, 513, 398

Shaw Greening, L., Barnard, R., Kolb, U., Tonkin, C., & Osborne, J. P. 2009, A&A, 495, 733

Siegel, M. H., Dotter, A., Majewski, S. R., et al. 2007, ApJL, 667, L57

Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton, NJ: Princeton Univ. Press)

Spitzer, L., Jr., & Hart, M. H. 1971, ApJ, 164, 399

Strader, J., Caldwell, N., & Seth, A. C. 2011, AJ, 142, 8

Supper, R., Hasinger, G., Lewin, W. H. G., et al. 2001, A&A, 373, 63

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

Trudolyubov, S., Kotov, O., Priedhorsky, W., Cordova, F., & Mason, K. 2005, ApJ, 634, 314

Trudolyubov, S., & Priedhorsky, W. 2004, ApJ, 616, 821

Trudolyubov, S. P., & Priedhorsky, W. C. 2008, ApJ, 676, 1218

van der Klis, M. 1998, AdSpR, 22, 925

Verbunt, F. 2003, in ASP Conf. Ser. 296, New Horizons in Globular Cluster Astronomy, ed. G. Piotto et al. (San Francisco, CA: ASP), 245

Verbunt, F., & Freire, P. C. 2014, A&A, 561, A11

Verbunt, F., & Hut, P. 1987, in IAU Symp. 125, The Origin and Evolution of Neutron Stars, ed. D. J. Helfand & J.-H. Huang (Dordrecht: Reidel), 187

Verbunt, F., Wijers, R. A. M. J., & Burm, H. M. G. 1990, A&A, 234, 195

White, N. E., & Angelini, L. 2001, ApJL, 561, L101

Wijnands, R., & van der Klis, M. 1998, Natur, 394, 344

Yungelson, L. R., Nelemans, G., & van den Heuvel, E. P. J. 2002, A&A, 388, 546