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

# Anti-glitches in the Ultraluminous Accreting Pulsar NGC 300 ULX-1
Observed with NICER

Paul S. Ray<fnref n="1" />, Sebastien Guillot<fnref n="2" />,<fnref n="3" />, Wynn C. G. Ho<fnref n="4" />,<fnref n="5" />, Matthew Kerr<fnref n="1" />, Teruaki Enoto<fnref n="6" />, Keith C. Gendreau<fnref n="7" />,
Zaven Arzoumanian<fnref n="7" />, Diego Altamirano<fnref n="8" />, Slavko Bogdanov<fnref n="9" />, Robert Campion<fnref n="10" />, Deepto Chakrabarty<fnref n="11" />,
Julia S. Deneva<fnref n="12" />, Gaurava K. Jaisawal<fnref n="13" />, Robert Kozon¹⁴, Christian Malacaria<fnref n="15" />,<fnref n="16" />, Tod E. Strohmayer<fnref n="17" />, and
Michael T. Wolff<fnref n="1" />

## Abstract

We present evidence for three spin-down glitches (or “anti-glitches”) in the ultraluminous accreting X-ray pulsar NGC 300 ULX-1, in timing observations made with the Neutron Star Interior Composition Explorer. Our timing analysis reveals three sudden spin-down events of magnitudes $\Delta\nu = -23, -30$, and $-43$ µHz (fractional amplitudes $\Delta\nu/\nu = -4.4, -5.5$, and $-7.7 \times 10^{-4}$). We determined fully phase-coherent timing solutions through the first two glitches, giving us high confidence in their detection, while the third candidate glitch is somewhat less secure. These are larger in magnitude (and opposite in sign) than any known radio pulsar glitch. This may be caused by the prolonged rapid spin up of the pulsar, causing a sudden transfer of angular momentum between the superfluid and non-superfluid components of the star. We find no evidence for profile or spectral changes at the epochs of the glitches, supporting the conclusion that these are due to the same process as in normal pulsar glitches, but in reverse.

**Key words:** pulsars: general – pulsars: individual (NGC 300 ULX-1) – X-rays: binaries

## 1. Introduction

Ultraluminous X-ray sources (ULXs) are extragalactic, bright X-ray point sources with super-Eddington luminosities for isotropic emission due to accretion onto a neutron star or stellar-mass black hole (see Kaaret et al. 2017, for a recent review). Following the discovery of X-ray pulsations from a handful of ULXs, with periods ranging from 0.42 s to ~31 s (Bachetti et al. 2014; Fürst et al. 2016; Israel et al. 2017a, 2017b), it has been suggested that the central engines of most ULXs could be highly magnetized accreting neutron stars, based on theoretical (King & Lasota 2016; Mushtukov et al. 2017) or observational arguments (Koliopanos et al. 2017; Walton et al. 2018b). Further evidence comes from the claimed discovery of cyclotron resonant scattering features (CRSFs) in two ULXs (Brightman et al. 2018; Walton et al. 2018a), albeit one of them with no detected pulsations.

The ULX in the nearby spiral galaxy NGC 300 ($D = 1.87 \pm 0.12$ Mpc; Rizzi et al. 2006) was mistakenly classified as a supernova (SN 2010da) upon its discovery (Monard 2010). Follow up Chandra observations showed that the system

is instead a high-mass X-ray binary undergoing an outburst (Binder et al. 2011). Using observations obtained in 2016 December with XMM-Newton and NuSTAR, Carpano et al. (2018b) discovered X-ray pulsations at a frequency of $\nu = 31.6$ mHz ($P = 31.6$ s, first reported in Carpano et al. 2018a), thus revealing the compact object as a neutron star. They also found that the pulsar is spinning up very rapidly with $\dot{\nu} = 5.57 \times 10^{-10}$ Hz s$^{-1}$.

This article reports on the regular monitoring of NGC 300 ULX-1 with the Neutron Star Interior Composition Explorer (NICER; Gendreau et al. 2016), and the discovery of three apparent spin-down glitches that occurred over a span of about 120 days. In the many hundreds of radio, X-ray, and γ-ray pulsars monitored for glitches, all of the hundreds of detected glitches were spin-up glitches ($\Delta\nu > 0$; see, e.g., Espinoza et al. 2011; Yu et al. 2013; Fuentes et al. 2017). Among magnetars, observations reveal approximately twenty glitches, consisting of both spin-up and spin-down ones (see, e.g., Dib & Kaspi 2014). The largest magnetar spin-down glitch (with $\Delta\nu = -21$ µHz) is inferred to possibly have occurred during a giant flare of SGR 1900+14 (Woods et al. 1999), while the second largest was over a hundred times smaller

<footnote n="1">¹Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375, USA; paul.ray@nrl.navy.mil</footnote>

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

<footnote n="3">³Université de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA</footnote>

<footnote n="5">⁵Mathematical Sciences, Physics and Astronomy, and STAG Research Centre, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="6">⁶The Hakubi Center for Advanced Research and Department of Astronomy, Kyoto University, Kyoto 606-8302, Japan</footnote>

<footnote n="7">⁷X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="8">⁸Physics and Astronomy, University of Southampton, Southampton, Hampshire SO17 1BJ, UK</footnote>

<footnote n="9">⁹Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="10">¹⁰Chesapeake Aerospace, Grasonville, MD 21638, USA</footnote>

<footnote n="11">¹¹MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="12">¹²George Mason University, resident at the Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="13">¹³National Space Institute, Technical University of Denmark, Elektrovej 327-328, DK-2800 Lyngby, Denmark
¹⁴Code 681, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="15">¹⁵NASA Marshall Space Flight Center, NSSTC, 320 Sparkman Drive, Huntsville, AL 35805, USA¹⁸</footnote>

<footnote n="16">¹⁶Universities Space Research Association, NSSTC, 320 Sparkman Drive, Huntsville, AL 35805, USA</footnote>

<footnote n="17">¹⁷Astrophysics Science Division and Joint Space-Science Institute, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

*Received 2018 December 7; revised 2019 May 23; accepted 2019 May 26; published 2019 July 16*</footnote>

<footnote n="18">¹⁸NASA Postdoctoral Fellow.</footnote>
(Δν ~ -0.09 μHz) and accompanied by radiative changes (Archibald et al. 2013). The glitches described above all occur in nonaccreting pulsars, whose long-term spin frequencies decrease with time (ν̇ < 0) due to rotational energy loss from electromagnetic dipole radiation (Pacini 1968; Gunn & Ostriker 1969). For pulsars that accrete matter from a binary companion, the torque due to accretion can cause the neutron star to spin up or spin down and is usually much stronger than the dipole radiation torque. There are only a few reported glitches in accreting pulsars; two spin-up glitches with Δν ≈ 1 μHz, each occurring during an accretion outburst episode (Galloway et al. 2004; Serim et al. 2017), and possibly two spin-down glitches (Klochkov et al. 2009).

## 2. Observations

In response to a series of Astronomer's Telegrams reporting the discovery of pulsations and the period evolution of NGC 300 ULX-1 (Bachetti et al. 2018a, 2018b; Carpano et al. 2018a; Grebenev & Mereminskiy 2018a, 2018b; Kennea 2018; Vasilopoulos et al. 2018a), NICER observed NGC 300 ULX-1 on 2018 February 6 for 1.6 ks. Following that, the source entered NICER's solar exclusion zone and became unobservable until 2018 May 1, when NGC 300 ULX-1 was again sufficiently far from the Sun. NICER then started a monitoring campaign on the source, primarily designed to discover evidence of an orbital period. NICER generally has a visibility window of a few hundred to ~2000 s in one segment of an orbit, and observations were generally made during 1–10 consecutive orbital segments per day (most often 2, 3, or 4). Here we present an analysis of NICER ObsIDs 1034200101 (2018 February 6) through 1034200198 (2018 October 20).

The X-ray timing instrument of NICER provides ≈1900 cm² of collecting area at 1.5 keV thanks to its 56 coaligned X-ray concentrating optics paired with single-pixel silicon drift detectors (SDDs) working in the 0.2–12 keV energy band (52 detectors are operational on orbit). The GPS time tagging of photons permits obtaining an accuracy better than 100 ns root mean square (rms; LaMarr et al. 2016; Prigozhin et al. 2016).

The data were processed following standard procedures using the NICER data analysis software NICERDAS version 4, together with HEASOFT version 6.24. The standard filtering criteria define the good time intervals for which the NICER pointing offset is <0°:015 from the nominal source position; the source is at least 30° from the Earth limb (at least 40° in the case of a Sun-illuminated Earth), and intervals when NICER was not within the boundaries of the South Atlantic Anomaly. In a few observations, we masked “hot” detectors exhibiting out-of-family count rates. Also, for ObsID 1034200101, we had to remove MPU3¹⁹ because of extreme count rates caused by solar loading on the detectors, and relaxed the filtering for the bright Earth to 30°. For each observation, we selected events in the 0.4–12 keV band and barycentered the data using the NASA/JPL DE405 solar system ephemeris and the coordinates of the ULX from Chandra imaging data (R.A. = 13°:770208, decl. = -37°:695417; Binder et al. 2011). All barycentric times and frequencies in this work use the Barycentric Dynamical Time (TDB) timescale. Some of the observations exhibited background flares when the International Space Station (ISS) was at extreme northern or southern

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Phaseogram of NGC 300 ULX-1 in the 0.4–12 keV band, from the NICER observations around MJD 58246.</figcaption></figure>

latitudes. We applied a count rate cut of 7.5 c s⁻¹ on non-overlapping 8 s intervals to exclude the times of such flares.

## 3. Timing

The timing analysis demonstrates the power of NICER for timing studies of faint X-ray pulsars. NICER's flexible scheduling enabled monitoring of the source in a way that maximized the ability to do precise timing. We scheduled observations to cover two or more consecutive ISS orbits, spreading a relatively small amount of observing time over a duration long enough to make a precise frequency measurement, but not so long that the window function caused significant ambiguity in the pulse counting between data segments. These precise frequency measurements allowed us to bootstrap a fully coherent timing solution. In the following, we present the timing analysis of the pulsations, working from coarse frequency measurements to phase-coherent timing.

For the initial timing analysis, we binned the data at 1/128 s and used the PRESTO software package (Ransom et al. 2002) to perform a periodicity search over a range of possible accelerations. In essentially every observation until 2018 August 22, when the source seems to have faded significantly, the ULX pulsations (at around 53–57 mHz) were detected as the highest significance peak coming out of the search. This allowed us to confirm that the period evolution was continuing as previously reported. In addition, we confirmed the broad Gaussian pulse profile of this ULX reported previously (Carpano et al. 2018b). Figure 1 shows the pulse profile from

<footnote n="19">¹⁹ MPU: Measurement Power Unit; each MPU connects 8 of the focal plane modules, each containing one SDD.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Pulse frequency evolution of NGC 300 ULX-1 as measured with NICER from each short, nearly contiguous data segment. The line is the fitted model as described in the text, and the bottom panel shows the residuals to the fit. Note that glitches are not readily apparent here because of the poor precision of the individual period measurements.</figcaption></figure>

MJD 58246, NICER's densest day of observations. The profile is reasonably well modeled as a Gaussian with a full width at half maximum (FWHM) of 0.35 cycles.

To go from initial periodicity detections to precise period measurements, we employed a maximum likelihood technique to measure pulsed frequencies and significance. In this case, given a selection of photons from a segment of the observation, we searched over a grid of trial $ν$ values. For each trial we computed the pulse phase of each photon according to a model with the trial $ν$ (and $\dot{ν}$ fixed at the long-term average, $4.16 \times 10^{-10}$ Hz s$^{-1}$). We modeled the probability distribution of this set of phases using a Gaussian template with a FWHM 0.35 of the pulse period and subsequently maximized the likelihood to determine the optimal values of the Gaussian peak position and amplitude. By profiling over these nuisance parameters, we arrived at a profile likelihood for $ν$; see Kerr et al. (2015) and Ray et al. (2011) for a discussion of similar likelihood techniques. This finally yielded a maximum likelihood estimate for the pulse frequency at that observation, as well as an estimate of the measurement uncertainty determined by the frequency range over which the log likelihood decreased by 0.5.

The rapid and noisy period evolution, combined with the window function of the observations, made precision timing difficult, so we developed a bootstrapping technique to make precise pulse frequency measurements. The first step was making likelihood frequency measurements for each nearly contiguous segment of NICER data (i.e., with no gaps greater than 1000 s and a total span less than 3000 s). These measurements are low-precision, but unambiguous. We note that the last significant pulsed detection was at MJD 58352.5 (2018 August 22; ObsID 1034200183). All subsequent observations show no strong evidence of pulsations from NGC 300 ULX-1; this includes an 800 s observation at MJD 58351.5, followed by a gap of a month, then 21 observations between MJD 58386 and 58411.<fnref n="20" />

The measured frequencies for each segment of NICER data are shown in Figure 2. Fitting the pulse frequency evolution to a simple model with two frequency derivatives gives the best fit $ν = ν_0 + \dot{ν}(t - t_0) + \frac{1}{2}\dot{ν}^2(t - t_0)^2$, where $ν_0 = 53.350(5)$ mHz, $\dot{ν} = 4.16(5) \times 10^{-10}$ Hz s$^{-1}$, $\dot{ν} = -2.5(12) \times 10^{-18}$ Hz s$^{-2}$, and $t_0 = \text{MJD}(\text{TDB}) 58243.515$. Here, and elsewhere in the paper, the numbers in parentheses are the $1\sigma$ uncertainty in the last digit. If the spin down proceeds as $\Omega \propto \Omega^n$ (where $\Omega$ is the angular velocity), the “braking index” $n$ is defined as $n = \nu\dot{\nu}/\dot{\nu}^2$. For NGC 300 ULX-1, $n$ is $-0.8(4)$. The small magnitude of $n$ indicates that we are likely seeing the secular spin up of the system and the $\dot{ν}$ is not strongly contaminated by timing noise or changing Doppler shift from an orbit. An analysis of the long-term spin up in terms of accretion torque theory is presented in Vasilopoulos et al. (2018b). The fitting procedure also yields a pulsed count rate for each observation; these are shown in Figure 3.

To make more precise frequency measurements, we repeated the likelihood analysis above, with expanded observation segments with a maximum span of 15,000 s (limiting the duration to three NICER orbits). In this stage, there are potential pulse counting ambiguities due to the long gaps between orbits in some observations. However, in nearly all cases, measurements that suffered from an off-by-one counting error between observation subsets were very clearly identified from the short-segment likelihood frequency estimates and the known position of the sidelobes from the observation window function.

These precise frequency measurements, excluding the single observation in 2018 February, are shown in Figure 4, where a constant frequency derivative has been subtracted to allow the fine structure to be seen. This reveals relatively smooth frequency evolution that shows some evidence for several negative frequency steps, which we identify as candidate glitches. However, the sampling and large uncertainties in some of the measurements, coupled with the rapid underlying spin evolution, makes this result tentative, requiring confirmation from the phase-coherent analysis presented in Section 4.

We attempted to fit the frequency evolution with a secular evolution up through $\dot{ν}$ and three glitches. The residuals to this fit are shown in Figure 5. There are clearly some unmodeled residuals, particularly in the high-precision measurements between glitches 1 and 2, and glitches 2 and 3. This could

<footnote n="20">²⁰ While this paper was in review, monitoring continued and the source rebrightened allowing some additional pulsation detections after MJD 58450. These were not sufficient to establish a coherent timing solution or to search for additional glitches. For an analysis of those data and the implications of the spin evolution, see Vasilopoulos et al. (2019).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Pulsed count rate for each NICER observation where pulsations from NGC 300 ULX-1 were significantly detected. The count rates near MJD 58155 have been scaled by 7/6 to account for the exclusion of MPU3 in the analysis. Observations with at least 500 s of exposure where no significant pulsation was detected are marked with red × symbols.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** NICER pulse frequency measurements for NGC 300 ULX-1, after subtraction of a simple model with a constant frequency derivative. The candidate glitch epochs are marked with red dashed lines.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Residuals to the best-fit three-glitch model for the NGC 300 ULX-1 frequency evolution. The upper (orange) points are the same data (offset by 2.0 × 10⁻⁴ Hz) with the glitch magnitudes set to 0, to show the location and amplitude of the glitches. For reference, the time ranges segments A–F used in the coherent analysis in Section 4 are also shown.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Pulse profiles (0.4–12 keV) from six different sections of data. Each profile has the mean subtracted and is normalized by the standard deviation. Two cycles are shown. The red points with error bars on the second cycle are the residuals of each profile to a model profile fitted to all of the data.</figcaption></figure>

<table><caption>**Table 1**
Incoherent Timing Parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Epoch (MJD)</td><td>58243.515</td></tr><tr><td>ν (mHz)</td><td>53.3533(7)</td></tr><tr><td>ν̇ (Hz s<sup>-1</sup>)</td><td>4.40(1) × 10<sup>-10</sup></td></tr><tr><td>ν̇̇ (Hz s<sup>-2</sup>)</td><td>-3.4(9) × 10<sup>-18</sup></td></tr><tr><td>ν<sup>(3)</sup> (Hz s<sup>-3</sup>)</td><td>-1.1(3) × 10<sup>-24</sup></td></tr><tr><td>Glitch 1 Epoch (MJD)</td><td>58243.515</td></tr><tr><td>Glitch 1 Δν (Hz)</td><td>-2.1(1) × 10<sup>-5</sup></td></tr><tr><td>Glitch 2 Epoch (MJD)</td><td>58265.5</td></tr><tr><td>Glitch 2 Δν (Hz)</td><td>-3.8(2) × 10<sup>-5</sup></td></tr><tr><td>Glitch 3 Epoch (MJD)</td><td>58333.5</td></tr><tr><td>Glitch 3 Δν (Hz)</td><td>-3.7(5) × 10<sup>-5</sup></td></tr></tbody></table>

potentially be modeled as frequency derivative steps at the glitch, or exponentially decaying glitch “recoveries.” Both of these phenomena are observed in rotation-powered pulsar glitches, but it is not clear which applies here, so we keep the model as simple as possible. The model parameters are shown in Table 1. We note that these model parameters yield a braking index $n = -0.9(2)$.

### 3.1. *Search for Profile Variations*

In most radio pulsar glitches the pulse profile does not change at the epoch of the glitch, though significant profile changes are observed at the epoch of some magnetar glitches (e.g., see Archibald et al. 2013) and radiative changes were

seen after a glitch in one high magnetic field pulsar PSR J1119 −6127 (Weltevrede et al. 2011). To help constrain the models for the NGC 300 ULX-1 glitches, we looked for profile variations. To do this in a timing-model-independent way, we considered observations where the pulsations were strongly detected (6σ or more). For each segment, we fitted a single Gaussian template to determine phase 0 and generated a set of phases using that as a reference. We then collected all of these phases from six different sections of the NICER monitoring observation, splitting the data at the glitches, and at the first significant data gap after each glitch. We binned these phases into 40-bin histograms, subtracted the mean and divided by the standard deviation, to account for the changing background and evolving pulsed count rate. These six pulse profiles are displayed in Figure 6. We also fitted a three-Gaussian template model to the complete data set and displayed the residuals of each profile to that mean shape. We saw no significant pulse profile changes during our observations.

### 4. Coherent Timing Analysis

Using the observed frequency evolution, we computed pulse times of arrival (TOAs) for each 150 s segment of NICER data (300 s for segment F, as defined below, where the count rate was significantly lower), using the unbinned maximum likelihood technique described in Ray et al. (2011). The reference phase for each TOA is the peak of the Gaussian pulse. Typical TOA uncertainties are 0.5 s, or about 3% of a pulse period.
<table><caption>**Table 2**
Parameters for Glitch 1 (Segment A)</caption><thead><tr><th colspan="2">Fit and Data Set</th></tr></thead><tbody><tr><td>MJD range</td><td>58238.9–58249.9</td></tr><tr><td>Data span (d)</td><td>10</td></tr><tr><td>Number of TOAs</td><td>385</td></tr><tr><td>rms timing residual (s)</td><td>0.607</td></tr><tr><td>Reduced χ<sup>2</sup> value</td><td>1.6</td></tr><tr><th colspan="2">Fixed Quantities</th></tr><tr><td>R.A., α (hh:mm:ss, J2000)</td><td>00:55:04.86</td></tr><tr><td>decl., δ (dd:mm:ss, J2000)</td><td>-37:41:43.7</td></tr><tr><td>Frequency second derivative, ṷ (Hz s<sup>-2</sup>)</td><td>0</td></tr><tr><td>Epoch of frequency determination (MJD)</td><td>58243.5</td></tr><tr><td>Glitch epoch (MJD)</td><td>58243.509876</td></tr><tr><th colspan="2">Measured Quantities</th></tr><tr><td>Pulse frequency, ν (mHz)</td><td>53.35285(19)</td></tr><tr><td>Frequency first derivative, ṷ (Hz s<sup>-1</sup>)</td><td>4.358(13) × 10<sup>-10</sup></td></tr><tr><td>Glitch Δν</td><td>-2.340(20) × 10<sup>-5</sup></td></tr><tr><td>Glitch Δν̇</td><td>1.41(13) × 10<sup>-11</sup></td></tr><tr><th colspan="2">Derived Quantities</th></tr><tr><td>Δν/ν̇</td><td>-4.39(4) × 10<sup>-4</sup></td></tr><tr><td>Δν̇/ν̇</td><td>3.2(3) × 10<sup>-2</sup></td></tr></tbody></table>

Because of the significant red noise, the rapid spin evolution, and the gaps between observations, generating a phase-connected timing model for the full data set is not possible. However, we were able to make phase-connected timing models for several contiguous sections of data (which we label as segments A–F), including spans that cover the first (segment A) and second (segment B) glitches. For later segments (C–F) of observations, we fitted simple models with $ν$ and $ν̇$. We use TEMPO2 (Hobbs et al. 2006) to fit timing models to the TOAs from each segment. The timing model parameters for each segment are presented in Tables 2–4, and the modeled spin-frequency evolution is displayed in Figure 7. The TOA residuals for each model are presented in Figures 8–10.

We considered the possibility of a glitch between segments C and D, but rejected that hypothesis since the extrapolated fits do not appear to require one. For segments E and F, we were able to get a good measurement of the frequency step at the glitch.

The coherent timing models shown in Figures 8–10 clearly demonstrate that there are sudden pulse frequency changes. However, it is conceivable that such features could be just a manifestation of the red noise in the system. To test this, we simulated red noise residuals. Most accretion-powered pulsars show white torque noise (Bildsten et al. 1997), which gives a pulse phase power spectrum with an index of −4, while some are even steeper. We made multiple realizations of residuals with red noise spectral indices of −3, −4, and −5, with normalizations set to match the rms of the observed residuals. The simulations with index −3 sometimes showed cusp-like features, but never with the large difference in scale between the red noise wiggles and the glitch residuals that we see in Figures 8 and 10. For more realistic red noise spectral indices of −4 and −5, the red noise is much smoother and extremely unlikely to be mistaken for a glitch. We therefore conclude that the timing features we observe are best characterized as glitches.

<table><caption>**Table 3**
Parameters for Glitch 2 (Segment B)</caption><thead><tr><th colspan="2">Fit and Data Set</th></tr></thead><tbody><tr><td>MJD range</td><td>58252.5–58278.7</td></tr><tr><td>Data span (d)</td><td>26</td></tr><tr><td>Number of TOAs</td><td>235</td></tr><tr><td>rms timing residual (μs)</td><td>1.36</td></tr><tr><td>Reduced χ<sup>2</sup> value</td><td>7.0</td></tr><tr><th colspan="2">Fixed Quantities</th></tr><tr><td>R.A., α (hh:mm:ss, J2000)</td><td>00:55:04.86</td></tr><tr><td>decl., δ (dd:mm:ss, J2000)</td><td>-37:41:43.7</td></tr><tr><td>Frequency second derivative, ṷ (Hz s<sup>-2</sup>)</td><td>0</td></tr><tr><td>Epoch of frequency determination (MJD)</td><td>58265</td></tr><tr><td>Glitch epoch (MJD)</td><td>58265.25263</td></tr><tr><th colspan="2">Measured Quantities</th></tr><tr><td>Pulse frequency, ν (mHz)</td><td>54.12969(9)</td></tr><tr><td>First derivative of pulse frequency, ṷ (s<sup>-2</sup>)</td><td>4.1603(17) × 10<sup>-10</sup></td></tr><tr><td>Glitch Δν</td><td>-2.990(14) × 10<sup>-5</sup></td></tr><tr><td>Glitch Δν̇</td><td>1.24(3) × 10<sup>-11</sup></td></tr><tr><th colspan="2">Derived Quantities</th></tr><tr><td>Δν/ν̇</td><td>-5.52(2) × 10<sup>-4</sup></td></tr><tr><td>Δν̇/ν̇</td><td>2.98(6) × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Piecewise coherent spin-frequency model for NGC 300 ULX-1. Each color represents a segment of data with a coherent timing model fit (denoted by the letters).</figcaption></figure>

## 5. Spectral Characteristics

Spectral analysis was performed with XSPEC version 12.10 (Arnaud 1996), and with the NICER response files version 0.06. For these analyses, the data filtering is more stringent than that of the timing analyses presented above. This is in an effort to minimize the non-astrophysical background that affects the spectra. First, we limit the energy range to 0.4–10 keV because of significant optical loading affecting channels below 0.4 keV, especially when NICER is observing at low Sun angles. In addition to the standard filtering described in Section 2, we also exclude the detectors with DET_ID 14, 34, and 54, which often suffer from optical loading. Finally, we exclude periods with the conditions: COR_SAX < 1.5, FPM_OVERONLY_COUNT > 1.0, and FPM_OVERONLY_COUNT > (1.52 × COR_SAX<sup>-0.633</sup>). This last condition is an empirical relation obtained from a large number of “blank sky” pointings to characterize NICER background over a wide range of orbital environments. Together, these remove time intervals with high particle background.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Coherent timing of glitch 2 (segment B). Left: residuals to a timing model including the glitch. Substantial red noise is visible in the residuals. Right: residuals to the same model with $\Delta\nu$, $\Delta\dot{\nu}$ set to 0 to show the effect of the glitch.</figcaption></figure>

<table><caption>**Table 4**

Coherent Timing Model Parameters for Segments C–F</caption><thead><tr><th colspan="5">Fit and Data Set</th></tr><tr><th>Sequence</th><th>C</th><th>D</th><th>E</th><th>F</th></tr></thead><tbody><tr><td>MJD range</td><td>58292.6–58296.0</td><td>58310.4–58318.4</td><td>58327.2–58332.3</td><td>58334.2–58353.2</td></tr><tr><td>Data span (d)</td><td>3.4</td><td>8</td><td>5.1</td><td>19</td></tr><tr><td>Number of TOAs</td><td>28</td><td>89</td><td>39</td><td>76</td></tr><tr><td>rms resid. (s)</td><td>0.93</td><td>0.88</td><td>0.73</td><td>2.0</td></tr><tr><td>Reduced χ<sup>2</sup></td><td>3.0</td><td>2.3</td><td>1.5</td><td>15.5</td></tr><tr><th colspan="5">Fixed Quantities</th></tr><tr><td>R.A., α</td><td></td><td></td><td>00:55:04.86</td><td></td></tr><tr><td>decl., δ</td><td></td><td></td><td>-37:41:43.7</td><td></td></tr><tr><td>ν̇ (s<sup>-3</sup>)</td><td></td><td></td><td>0</td><td></td></tr><tr><td>Epoch (MJD)</td><td>58298</td><td>58298</td><td>58333</td><td>58333</td></tr><tr><th colspan="5">Measured Quantities</th></tr><tr><td>ν (mHz)</td><td>55.3280(8)</td><td>55.3295(5)</td><td>56.5046(3)</td><td>56.46146(13)</td></tr><tr><td>ν̇ (Hz s<sup>-1</sup>)</td><td>4.25(3) × 10<sup>-10</sup></td><td>3.877(3) × 10<sup>-10</sup></td><td>3.838(9) × 10<sup>-10</sup></td><td>3.8299(14) × 10<sup>-10</sup></td></tr><tr><th colspan="5">Derived Quantities</th></tr><tr><td>Δν (Hz)</td><td>(D–C) 1.5(9) × 10<sup>-6</sup></td><td></td><td>(F–E) -4.32(3) × 10<sup>-5</sup></td><td></td></tr><tr><td>Δν̇/ν̇</td><td>2.7(16) × 10<sup>-5</sup></td><td></td><td>-7.66(6) × 10<sup>-4</sup></td><td></td></tr></tbody></table>

We generated a background spectrum specific to small time intervals using a library of spectra populated with NICER observations of seven blank sky regions, corresponding to the RXTE background regions (Jahoda et al. 2006). The background spectra we obtain are a combination of these blank sky spectra weighted according to space-weather observing conditions common to both the NGC 300 ULX-1 and background-field observations.

To each ObsID spectrum with its corresponding background, we fit a simple absorbed power law to study the spectral evolution over the several weeks of NICER observations of NGC 300 ULX-1. No obvious spectral changes can be determined at the epochs of the glitches. However, large variations of the background count rates estimated for each ObsID may hamper the precise study of the evolution.

Instead, we focus our study on the spectral evolution of the pulsed component only. For each ObsID, we extract spectra from the on-pulse phase range (0.75–1.25) and off-pulse phase range (0.25–0.75; see Figure 6). The off-pulse spectra, used as backgrounds for the on-pulse spectra, therefore permit subtraction of the non-astrophysical background and of the other X-ray sources within the NICER field of view. The typical flux of the pulsed component is $f_X \sim 2 \times 10^{-12} \text{ erg s}^{-1} \text{ cm}^{-2}$ (0.4–10 keV).

We fit each of these on-minus-off spectra with a simple absorbed power law. The best-fit parameters are shown in Figure 11, together with the reduced $\chi^2$ values. No significant spectral change is observed at the time of glitches 1 and 2. The observations around glitch 3 were shorter and provided low signal-to-noise spectra for a similar study. We therefore conclude that no detectable radiative changes were associated with the glitches.

## 6. Discussion

### 6.1. Characterization of the Glitches

We measured three glitches interrupting the rapid spin up of NGC 300 ULX-1, with magnitudes $\Delta\nu = -23, -30$, and $-43$ µHz. In the first two glitches, our best-fit coherent timing models include frequency derivative steps of $\Delta\dot{\nu} = 1.4$ and $1.2 \times 10^{-11} \text{ Hz s}^{-1}$, while in the third (more poorly sampled) glitch, we observe no significant $\Delta\dot{\nu}$.

With approximately daily sampling, we do not catch any of the glitches as they happen. We can generally only constrain the epoch of the glitch to about $\pm 0.5$ day accuracy, unless we impose the constraint that the pulse phase be continuous through the glitch, which gives a couple of more precise possible glitch epochs. There are significant unmodeled
<figure><img src="image_11.png" /><figcaption>Figure 9. Residuals to the coherent timing models for segments C, D, E, and F.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Coherent timing of glitch 1 (segment A). Left: residuals to a timing model including the glitch. Right: residuals to the same model with $\Delta\nu$, $\Delta\bar{\nu}$ set to 0 to show the effect of the glitch.</figcaption></figure>

residuals (red noise) evident in the coherent fits to the glitches 2 and 3; however, they are much smaller than the inferred glitch-induced component. We fit for a $\Delta\bar{\nu}$ at the glitch, but the nature of the red noise could include glitch recovery, accretion torque noise, and Doppler shifts from the orbit of the neutron star. We are unable to distinguish among these possibilities at this point. We note that the residuals look periodic in Figure 8. If fitted with an orbit, one gets a period of $P_B = 7.3$ day, with a semimajor axis ($a_1 \sin i)/c = 1.8$ lt-s. This does not seem like a plausible orbit for the neutron star as it implies a median companion mass of only $0.07 M_\odot$ (for an inclination of 60°), when the companion mass is known to be $20 M_\odot$ or larger. In addition, this oscillation does not fit the other data sets with the same period and phase, so we conclude that it is not due to an orbit.

## 6.2. Nature of the Glitches

The glitches measured in NGC 300 ULX-1 are unique among glitches observed to date. The lack of strong evidence

for any radiative change associated with the glitch in NGC 300 ULX-1 suggests a glitch mechanism internal to the neutron star, rather than a magnetospheric mechanism. Even in the case of magnetar glitches, where only some glitches are accompanied by radiative changes, they are thought to be caused by a source within the neutron star (Dib & Kaspi 2014). The internal mechanism that produces glitches in nonaccreting pulsars is thought to be due to a superfluid in the pulsar's inner crust and possibly core. The free neutrons that are present there are expected to be in a superfluid state because the typical crust temperature is well below the critical temperature for superfluidity. Unlike normal matter, superfluid matter rotates by forming vortices whose areal density determines the spin rate of the superfluid. These vortices are normally pinned to the lattice of crust nuclei, and as a result, while the rest of the star slows down owing to electromagnetic dipole radiation, the superfluid does not. A lag develops between the stellar spin rate and that of the superfluid. When this lag reaches a critical value, superfluid vortices unpin and transfer angular
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Pulsed component spectral evolution spanning segments A and B, over glitches 1 and 2. Each spectrum was fitted with an absorbed power law of index $\Gamma_{PL}$ and normalization $N_{PL}$. Vertical error bars are 90% confidence level uncertainties on the parameters, while the horizontal error bars represent the span of the ObsIDs. The dashed vertical lines show the epochs of glitches 1 and 2. The bottom panel shows the Kp space-weather index, an indicator of high backgrounds in NICER data.</figcaption></figure>

momentum to the rest of the star, causing the stellar spin rate to increase, i.e., a spin-up glitch (Anderson & Itoh 1975; see, e.g., Haskell & Melatos 2015, for a review; see left panel of Figure 12).

The above mechanism can still apply to the case of a spin-down glitch (see, e.g., Pines et al. 1980; see right panel of Figure 12), such as that seen here for NGC 300 ULX-1. Here, the pulsar spin frequency increases with time due to accretion (pre-glitch $\dot{v} = 4 \times 10^{-10}$ Hz s$^{-1}$), with a characteristic spin-up timescale $\tau_c \equiv \nu/2\dot{v} = 2$ yr. Because of this short timescale, we might expect a lag to build up quickly, where in our case most of the star spins faster than the superfluid. Once the lag is too great, the superfluid receives angular momentum from the rest of the star, and the star undergoes a spin-down glitch. This can be shown from the calculation of e.g., Andersson et al. (2012), by simply replacing the torque due to dipole radiation with that due to accretion. Such an interpretation provides a constraint on the superfluid fraction in the neutron star (Link et al. 1999; Andersson et al. 2012; Chamel 2013) and even potentially a determination of the neutron star mass (Ho et al. 2015, 2017). The former yields $I_{sf}/I \approx (m_n^{\text{eff}}/m_n)(-\sum_i \Delta\nu^i/\nu)(2\tau_c/t_{\text{obs}}) \approx 0.02(m_n^{\text{eff}}/m_n)(116 \text{ d}/t_{\text{obs}})$, where the superfluid and total moments of inertia are $I_{sf}$ and $I$, respectively, the average effective neutron mass

due to entrainment is $m_n^{\text{eff}}/m_n$ (~4.3; Chamel 2012), the summation is over each glitch, and the observation time is $t_{\text{obs}}$. The large glitch amplitudes ($|\Delta\nu| \approx 20-40$ µHz) and frequency derivative ($\dot{v} = 4 \times 10^{-10}$ Hz s$^{-1}$) of NGC 300 ULX-1 are comparable to those of the frequent large glitching pulsars Vela and PSR J0537−6910 and $\Delta\dot{v} \approx 10^{-11}$ Hz s$^{-1}$ are among the largest measured (Espinoza et al. 2011; Antonopoulou et al. 2018; Ferdman et al. 2018). This could suggest that NGC 300 ULX-1 has a similarly high glitch activity (Espinoza et al. 2011; Yu et al. 2013; Fuentes et al. 2017) if the underlying glitch mechanism is the same.

The rise time of a glitch is uncertain, due to unknown properties of superfluidity, but is thought to be >8 *P* (Sidery et al. 2010; Sourie et al. 2017). Such a limit is in accord with glitches seen in the Vela pulsar (Palfreyman et al. 2018) and suggests a short rise time of >150 s for NGC 300 ULX-1, well within the ~12 hr of our nearest post-glitch TOA (see Section 3).

### 6.3. Glitch Recurrence Timescale

We find three glitches in 116 days of monitoring, giving a rate of one glitch per 39 days. Alternatively, the two measured interglitch times are 22 and 68 days. Ducci et al. (2015) consider a model of the evolution of superfluid vortices that can produce a spin-down glitch. For this system, depending on the vortex pinning force, their model predicts a maximum amplitude of 2–20 µHz, glitch recovery timescale of $10^2-10^3$ s, and recurrence time of 40 days, although the observed spin-down $\Delta\dot{v}/\dot{v}$ seem to be much smaller than that estimated from this model, with the discrepancy possibly due to model extrapolations and uncertainties.

### 6.4. No Evidence for an Orbit

The initial motivation for these timing observations was to determine the orbit of the neutron star. To search for a short-period orbit, we used a Lomb–Scargle periodogram of the frequency residuals (Figure 5). No significant periodic Doppler shifts are detected.

If we assume the pulsar is a neutron star in orbit with a $30 M_{\odot}$ companion (suggested masses range from 20 to $40 M_{\odot}$; see Binder et al. 2016; Lau et al. 2016), for a typical inclination of $i = 60^{\circ}$, the mass function will be about $17.8 M_{\odot}$. From this, we can compute the expected amplitude of the orbital Doppler shifts, given a choice of orbital period, $P_B$. For a relatively short orbital period of a few days, the predicted amplitude would be easily detected in our data. For example, for $P_B = 7$ days, the amplitude of the sinusoidal Doppler modulation would be 0.05 mHz, which is much larger than the frequency residuals in Figure 5. Therefore, unless the system is very nearly face-on, short orbital periods with massive companions can be excluded. As described above, the residuals to the coherent fits may suggest a short orbit around a very low mass companion ($<0.1 M_{\odot}$), but we conclude that this is just a few cycles of a red noise process.

If, on the other hand, the orbital period is long (e.g., ~400 days, as suggested by Lau et al. 2016), the amplitude of the frequency modulation would be about 0.01 mHz, but on a longer timescale, so it would be absorbed in the $\dot{v}$ and $\ddot{v}$ fit to our frequency evolution. The expected $\dot{v}$ for a 400-day orbit with the same component masses is two orders of magnitude smaller than the observed value, so it would not contribute significantly. But, the $\ddot{v}$ magnitude would be up to $5 \times 10^{-19}$ Hz s$^{-2}$ (depending on
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left: schematic of spin-up glitch. Pulsar spin frequency (solid line) is observed to decrease due to energy loss from electromagnetic radiation, while a superfluid component within the star is pinned and does not spin down with the rest of the star (dashed line). When the superfluid component unpins and transfers angular momentum to the rest of the star, a spin-up glitch occurs. The pulsar then continues to spin down when the superfluid becomes pinned again. Right: schematic of spin-down glitch. Analogous to a spin-up glitch, except in this case the pulsar spin frequency (solid line) is observed to increase due to accretion torque and a spin-down glitch occurs when the superfluid gains angular momentum (dashed line).</figcaption></figure>

the orbital phase and inclination), which is a factor of 7 below the observed $\dot{v}$. Therefore, we conclude that a long orbital period could be contributing to the measured $\dot{v}$ but the majority of the observed frequency evolution is from accretion torques.

## 6.5. Magnetic Field and Energetics

Using an X-ray luminosity of $5 \times 10^{39}$ erg s$^{-1}$ (Carpano et al. 2018b) and standard accretion torque of $\dot{M}r_m^2\Omega_K(r_m)$, where $r_m$ is the approximate magnetosphere radius and $\Omega_K(r_m)$ is the Keplerian orbital angular frequency at $r_m$, we estimate that the magnetic field of NGC 300 ULX-1 is $\sim 1 \times 10^{12}$ G, which is in agreement with that determined from a possible detection of an electron CRSF at $E \sim 13$ keV (Walton et al. 2018a), although the field can be much higher if the pulsar is near spin equilibrium, where spin-up and spin-down torques are nearly equal. In the case of a proton CRSF, the inferred field would be $2 \times 10^{15}$ G/(1 + $z_g$). We note that the presence of this CRSF was called into question by Koliopanos et al. (2019), demonstrating that no spectral feature was required when modeling the XMM and NuSTAR spectra with a multicolor thermal emission model plus a nonthermal tail.

## 7. Summary

In this article, we presented evidence for two spin-down glitches (and a tentative third one) during the spin up of the ULX source NGC 300 ULX-1. Our NICER monitoring campaign allowed us, using phase-coherent timing analyses, to determine with high confidence that the magnitudes of these two glitches are the largest among known radio pulsar glitches. While the third glitch had an even larger magnitude ($\Delta\nu = -43$ µHz), its identification is somewhat less secure, as gaps in monitoring precluded precise determination of the glitch epoch and amplitude. We argue that these spin-down glitches are similar in nature to the spin-up glitches observed in spinning-down radio pulsars, i.e., the glitches are caused by sudden angular momentum transfer between a superfluid component and the rest of the neutron star. This is supported by the lack of detected radiative changes, which indicates a process internal to the neutron star. Finally, despite our dense

monitoring campaign, our timing solution did not permit determination of the orbit of the binary system, either because the system's inclination is near 0° or because the Doppler amplitude of a long period orbit (hundreds of days) is hidden in the frequency evolution ($\dot{v}$ and $\ddot{v}$) of the accretion torques in this system. Continued dense monitoring would permit studying the spin period evolution and accretion mechanisms of this intriguing source, as well as detecting more spin-down glitches to characterize the superfluid components inside the neutron star.

We thank the anonymous referee for useful comments and suggestions that improved this article. We thank Nishanth Anand (Thomas Jefferson High School for Science and Technology) for his work on the analysis for this paper during his SEAP summer internship at NRL. We thank Faith Huynh (MIT) for developing excellent NICER quicklook plotting software used in this work. We also thank Breanna Binder for helpful discussions. This work was supported by NASA through the NICER mission and the Astrophysics Explorers Program. NICER work at NRL is supported by NASA. S.G. is supported by the Centre National d'Etudes Spatiales (CNES). W.C.G.H. is partially supported by the Science and Technology Facilities Council (STFC) in the UK. This research has made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. This research has also made use of the NASA Astrophysics Data System (ADS) and the arXiv.

Facility: NICER.

Software: astropy (Astropy Collaboration et al. 2013), PINT (https://github.com/nanograv/pint), HEAsoft (https://heasarc.nasa.gov/lheasoft/), TEMPO2 (Hobbs et al. 2006), PRESTO (https://www.cv.nrao.edu/~sransom/presto/).

## ORCID iDs

Paul S. Ray https://orcid.org/0000-0002-5297-5278

Sebastien Guillot https://orcid.org/0000-0002-6449-106X

Wynn C. G. Ho https://orcid.org/0000-0002-6089-6836
Matthew Kerr https://orcid.org/0000-0002-0893-4073

Teruaki Enoto https://orcid.org/0000-0003-1244-3100

Diego Altamirano https://orcid.org/0000-0002-3422-0074

Slavko Bogdanov https://orcid.org/0000-0002-9870-2742

Deepto Chakrabarty https://orcid.org/0000-0001-8804-8946

Julia S. Deneva https://orcid.org/0000-0003-1226-0793

Gaurava K. Jaisawal https://orcid.org/0000-0002-6789-2723

Christian Malacaria https://orcid.org/0000-0002-0380-0041

Tod E. Strohmayer https://orcid.org/0000-0001-7681-5845

Michael T. Wolff https://orcid.org/0000-0002-4013-5650

## References

Anderson, P. W., & Itoh, N. 1975, Natur, 256, 25

Andersson, N., Glampedakis, K., Ho, W. C. G., & Espinoza, C. M. 2012, PhRvL, 109, 241103

Antonopoulou, D., Espinoza, C. M., Kuiper, L., & Andersson, N. 2018, MNRAS, 473, 1644

Archibald, R. F., Kaspi, V. M., Ng, C. Y., et al. 2013, Natur, 497, 591

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

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bachetti, M., Grefenstette, B. W., Walton, D. J., et al. 2018a, ATel, 11282, 1

Bachetti, M., Grefenstette, B. W., Walton, D. J., et al. 2018b, ATel, 11285, 1

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

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

Binder, B., Williams, B. F., Kong, A. K. H., et al. 2011, ApJL, 739, L51

Binder, B., Williams, B. F., Kong, A. K. H., et al. 2016, MNRAS, 457, 1636

Brightman, M., Harrison, F. A., Fürst, F., et al. 2018, NatAs, 2, 312

Carpano, S., Haberl, F., & Maitra, C. 2018a, ATel, 11158

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

Chamel, N. 2012, PhRvC, 85, 035801

Chamel, N. 2013, PhRvL, 110, 011101

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

Ducci, L., Pizzochero, P. M., Doroshenko, V., et al. 2015, A&A, 578, A52

Espinoza, C. M., Lyne, A. G., Stappers, B. W., & Kramer, M. 2011, MNRAS, 414, 1679

Ferdman, R. D., Archibald, R. F., Gourgouliatos, K. N., & Kaspi, V. M. 2018, ApJ, 852, 123

Fuentes, J. R., Espinoza, C. M., Reisenegger, A., et al. 2017, A&A, 608, A131

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

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

Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, Proc. SPIE, 9905, 99051H

Grebenev, S. A., & Mereminskiy, I. A. 2018a, ATel, 11174, 1

Grebenev, S. A., & Mereminskiy, I. A. 2018b, ATel, 11228, 1

Gunn, J. E., & Ostriker, J. P. 1969, Natur, 221, 454

Haskell, B., & Melatos, A. 2015, IJMPD, 24, 1530008

Ho, W. C. G., Espinoza, C. M., Antonopoulou, D., & Andersson, N. 2015, SciA, 1, e1500578

Ho, W. C. G., Espinoza, C. M., Antonopoulou, D., & Andersson, N. 2017, JPSCP, 14, 010805

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

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

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

Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, ApJS, 163, 401

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

Kennea, J. A. 2018, ATel, 11229, 1

Kerr, M., Ray, P. S., Johnston, S., Shannon, R. M., & Camilo, F. 2015, ApJ, 814, 128

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

Klochkov, D., Staubert, R., Postnov, K., Shakura, N., & Santangelo, A. 2009, A&A, 506, 1261

Koliopanos, F., Vasilopoulos, G., Buchner, J., Maitra, C., & Haberl, F. 2019, A&A, 621, A118

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

LaMarr, B., Prigozhin, G., Remillard, R., et al. 2016, Proc. SPIE, 9905, 99054W

Lau, R. M., Kasliwal, M. M., Bond, H. E., et al. 2016, ApJ, 830, 142

Link, B., Epstein, R. I., & Lattimer, J. M. 1999, PhRvL, 83, 3362

Monard, L. A. G. 2010, CBET, 2289, 1

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

Pacini, F. 1968, Natur, 219, 145

Palfreyman, J., Dickey, J. M., Hotan, A., Ellingsen, S., & van Straten, W. 2018, Natur, 556, 219

Pines, D., Shaham, J., Alpar, M. A., & Anderson, P. W. 1980, PThPS, 69, 376

Prigozhin, G., Gendreau, K., Doty, J. P., et al. 2016, Proc. SPIE, 9905, 99051I

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

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

Rizzi, L., Bresolin, F., Kudritzki, R.-P., Gieren, W., & Pietrzyński, G. 2006, ApJ, 638, 766

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

Sidery, T., Passamonti, A., & Andersson, N. 2010, MNRAS, 405, 1061

Sourie, A., Chamel, N., Novak, J., & Oertel, M. 2017, MNRAS, 464, 4641

Vasilopoulos, G., Haberl, F., Carpano, S., & Maitra, C. 2018a, ATel, 11179, 1

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

Vasilopoulos, G., Petropoulou, M., Koliopanos, F., et al. 2019, arXiv:1905.03740

Walton, D. J., Bachetti, M., Fürst, F., et al. 2018a, ApJL, 857, L3

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

Weltevrede, P., Johnston, S., & Espinoza, C. M. 2011, MNRAS, 411, 1917

Woods, P. M., Kouveliotou, C., van Paradijs, J., et al. 1999, ApJL, 524, L55

Yu, M., Manchester, R. N., Hobbs, G., et al. 2013, MNRAS, 429, 688