# The Discovery of the 528.6 Hz Accreting Millisecond X-Ray Pulsar MAXI J1816-195

Peter Bult, Diego Altamirano, Zaven Arzoumanian, Deepto Chakrabarty, Jérôme Chenevez, Elizabeth C. Ferrara, Keith C. Gendreau, Sebastien Guillot, Tolga Güver, Wataru Iwakiri, et al.

► To cite this version:

Peter Bult, Diego Altamirano, Zaven Arzoumanian, Deepto Chakrabarty, Jérôme Chenevez, et al. The Discovery of the 528.6 Hz Accreting Millisecond X-Ray Pulsar MAXI J1816-195. The Astrophysical journal letters, 2022, 935, 10.3847/2041-8213/ac87f9. insu-03863536

HAL Id: insu-03863536

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

Submitted on 21 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# The Discovery of the 528.6 Hz Accreting Millisecond X-Ray Pulsar MAXI J1816–195

Peter Bult<sup>1,2 ID</sup>, Diego Altamirano<sup>3 ID</sup>, Zaven Arzoumanian<sup>2</sup>, Deepto Chakrabarty<sup>4 ID</sup>, Jérôme Chenevez<sup>5 ID</sup>, Elizabeth C. Ferrara<sup>1,2 ID</sup>, Keith C. Gendreau<sup>2 ID</sup>, Sebastien Guillot<sup>6 ID</sup>, Tolga Güver<sup>7,8 ID</sup>, Wataru Iwakiri<sup>9 ID</sup>, Gaurava K. Jaisawal<sup>5 ID</sup>, Giulio C. Mancuso<sup>10,11 ID</sup>, Christian Malacaria<sup>12 ID</sup>, Mason Ng<sup>4 ID</sup>, Andrea Sanna<sup>13 ID</sup>, Tod E. Strohmayer<sup>14 ID</sup>, Zorawar Wadiasingh<sup>1,2,15 ID</sup>, and Michael T. Wolff<sup>16 ID</sup>

Received 2022 July 14; revised 2022 August 8; accepted 2022 August 8; published 2022 August 22

## Abstract

We present the discovery of 528.6 Hz pulsations in the new X-ray transient MAXI J1816–195. Using NICER, we observed the first recorded transient outburst from the neutron star low-mass X-ray binary MAXI J1816–195 over a period of 28 days. From a timing analysis of the 528.6 Hz pulsations, we find that the binary system is well described as a circular orbit with an orbital period of 4.8 hr and a projected semimajor axis of 0.26 lt-s for the pulsar, which constrains the mass of the donor star to 0.10–0.55 $M_\odot$. Additionally, we observed 15 thermonuclear X-ray bursts showing a gradual evolution in morphology over time, and a recurrence time as short as 1.4 hr. We did not detect evidence for photospheric radius expansion, placing an upper limit on the source distance of 8.6 kpc.

Unified Astronomy Thesaurus concepts: Low-mass x-ray binary stars (939); Millisecond pulsars (1062); X-ray bursts (1814)

## 1. Introduction

Accreting millisecond X-ray pulsars (AMXPs) are rapidly rotating neutron stars that accrete matter from a binary companion (see Patruno & Watts 2021; Di Salvo & Sanna 2022, for reviews). The characteristic that sets these systems apart from the wider population of low-mass X-ray binary systems is that they exhibit a coherent pulsation that directly tracks the millisecond stellar rotation period of the neutron star. Such pulsations are a useful diagnostic for the accreting neutron star. For instance, the precise waveform of the pulsations encodes information about the shape of the surface emission region and the neutron star compactness, and thus its equation of state (Poutanen & Gierlinski 2003), while tracking of the pulse arrival times allows for a high precision measurement of the neutron star spin and binary ephemeris, and may in principle be used to investigate the torques acting on these millisecond pulsars (Bildsten 1998; Psaltis & Chakrabarty 1999; see Burderi et al. 2006; Hartman et al. 2008; Patruno & Watts 2021; Di Salvo & Sanna 2022 for recent discussions of different torque mechanisms that may play a role).

An enduring challenge to the study of AMXPs is that they are relatively rare. Since the discovery of pulsations from SAX J1808.4–3658 in 1998 (Wijnands & van der Klis 1998), the population of AMXPs has grown at an average rate of under one per year. All the known AMXPs are X-ray transients, and the accretion-powered pulsations are only visible while the source is undergoing an X-ray outburst. These outbursts typically last only a few days to a few weeks, and are interspersed by several years to decades of inactivity (Lasota 2001; Hameury 2020). As such, the search for new accreting millisecond pulsars remains an important task.

In this Letter we present the discovery of 528.6 Hz pulsations from MAXI J1816–195 (henceforth MAXI J1816). This system was first discovered as a new X-ray transient with the MAXI Gas Slit Camera on 2022 June 7 (Negoro et al. 2022), and an initial source localization was provided by Swift shortly thereafter (Kennea et al. 2022a). Subsequent follow-up observations with the Neutron Star Interior Composition Explorer (NICER) on 2022 June 8 revealed the presence of 528.6 Hz pulsations (Bult et al. 2022a), identifying MAXI J1816 as an accreting millisecond pulsar.

Further monitoring with NICER revealed the 4.8 hr binary orbit of MAXI J1816 (Bult et al. 2022b). Additional Swift observations yielded an improved X-ray localization (Kennea et al. 2022b), after which likely counterparts to the X-ray source were identified in infrared (Kennea et al. 2022b), optical

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1"><sup>1</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="2"><sup>2</sup>Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="3"><sup>3</sup>Physics & Astronomy, University of Southampton, Southampton, Hampshire SO17 1BJ, UK</footnote>

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

<footnote n="5"><sup>5</sup>DTU Space, Technical University of Denmark, Elektrovej 327-328, DK-2800 Kongens Lynby, Denmark</footnote>

<footnote n="6"><sup>6</sup>Institut de Recherche en Astrophysique et Planétologie, UPS-OMP, CNRS, CNES, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="7"><sup>7</sup>Istanbul University, Science Faculty, Department of Astronomy and Space Sciences, Beyazıt, 34119, Istanbul, Turkey</footnote>

<footnote n="8"><sup>8</sup>Istanbul University Observatory Research and Application Center, Istanbul University 34119, Istanbul Turkey</footnote>

<footnote n="9"><sup>9</sup>Department of Physics, Faculty of Science and Engineering, Chuo University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo 112-8551, Japan</footnote>

<footnote n="10"><sup>10</sup>Instituto Argentino de Radioastronomía (CCT-La Plata, CONICET; CICPBA), C.C. No. 5, 1894 Villa Elisa, Argentina</footnote>

<footnote n="11"><sup>11</sup>Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque s/n, 1900 La Plata, Argentina</footnote>

<footnote n="12"><sup>12</sup>International Space Science Institute (ISSI), Hallerstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="13"><sup>13</sup>Dipartimento di Fisica, Università degli Studi di Cagliari, SP Monserrato-Sestu km 0.7, I-09042 Monserrato, Italy</footnote>

<footnote n="14"><sup>14</sup>Astrophysics Science Division and Joint Space-Science Institute, NASA’s Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="15"><sup>15</sup>Center for Research and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD 20771, USA</footnote>

<footnote n="16"><sup>16</sup>Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>
(de Martino et al. 2022), and radio (Beauchamp et al. 2022; Bright et al. 2022).

In the following we present the detailed analysis under-
pinning the NICER discovery of pulsations from MAXI J1816,
and the subsequent characterization of its orbit. We describe the
NICER monitoring campaign and provide a pulsar timing
analysis. Further, we present an analysis of fifteen thermo-
nuclear X-ray bursts observed from this system, and interpret
the various implications from both the bursts and pulsations on
the nature of MAXI J1816. An analysis of the (nonburst) X-ray
spectroscopy will be presented elsewhere.

## 2. Observations

We have monitored the 2022 outburst of MAXI J1816
extensively with NICER. Our observations began on 2022 June
7 and continued through 2022 July 5, at which time the source
intensity had decreased to the background level. These data are
available under the NICER ObsIDs starting with 520282 and
553301.

The source coordinates used for instrument pointing evolved
slightly over the course of the monitoring campaign. Our first
two observations were collected by pointing at the NICER
raster-scan coordinates reported by Bult et al. (2022a). The
subsequent monitoring program used the initial Swift localiza-
tion of Kennea et al. (2022a) until 2022 June 30, when we
switched to the corrected coordinates of Kennea et al. (2022b).
These coordinates are all <50" apart, which is much smaller
than the offset angle at which the response of NICER is
significantly affected.<fnref n="17" />

The data were processed using NICERDAS version 9 as
distributed with HEASOFT version 6.30. To account for the
small variations in the pointing coordinates, we filter the data
such that the angular offset relative to the pointing coordinates
is <54".<fnref n="18" /> Otherwise, we used standard filtering criteria; we
included only epochs that had an elevation angle >15°, a bright
Earth limb angle >30°, and were outside the South Atlantic
Anomaly (SAA).

After applying these filtering criteria, we were left with 92 ks
of clean exposure. We corrected the clean data to the solar
system barycenter using the DE-430 planetary ephemeris
(Folkner et al. 2014) and the source coordinates of Kennea
et al. (2022b).

Inspecting a 1 s time resolution light curve of both clean and
unfiltered data, we readily identify 15 thermonuclear (type I)
X-ray bursts. Ten of these were observed in full in the clean
data, while another three were observed during SAA passages.
The remaining two bursts were truncated by the boundaries of
the observations.

## 3. Results

We began our analysis by considering the evolution of the
outburst. We divided the data into segments of continuous
pointing, finding 93 such pointings across the 28 days of
monitoring with an exposure per pointing between
150 and 2250 s. For each pointing we calculated the average
count-rate in the 0.5–10 keV band, which we show in Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Outburst evolution of MAXI J1816 relative to 2022 June 7 (MJD 59737). The top panel shows the 0.5–10 keV count-rate as a function of time, with each point representing a single NICER pointing. The X-ray bursts have been removed from this data, but their onset times are marked with red crosses. The second panel shows the fractional pulse amplitude in the 0.5–10 keV range using triangles to indicate the 95% upper limits for nondetections. Finally, the two bottom panels give the phase residuals relative to a constant frequency and flux-bias adjusted timing model (see Table 1).</figcaption></figure>

3.1. Coherent Timing

We initially searched for the presence of pulsations by
dividing the first observation into 64 s segments and taking the
Fourier transform of each segment. Converting these trans-
forms to an averaged power spectrum revealed a high
significance (>6σ) pulse signal at 528.6 Hz.

To investigate the pulse signal in greater detail, we
determined the pulse frequency that optimized the $Z_1^2$ score
(Buccheri et al. 1983) for each pointing separately. The
resulting pulse frequencies were found to show a clear 4.8 hr
periodic oscillation in time, revealing the orbit of the binary
system. We fitted these frequency measurements using a
sinusoid to extract an initial estimate of the orbital period, $P_{\text{orb}}$,
the neutron star’s projected semimajor axis, $a_x \sin i$, and the
time of its passage through the ascending node, $T_{\text{asc}}$.

Based on the initial timing solution, we corrected the photon
arrival times for the Doppler delays of the binary motion and
then folded each pointing to a pulse profile. We fitted these
profiles using a constant plus two harmonically related
sinusoids, one fixed at the pulse frequency and the other at
twice that frequency. Either harmonic was deemed significant if
its measured amplitude was greater than three times the

<footnote n="17">¹⁷ https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/cal-recommend/</footnote>

<footnote n="18">¹⁸ The default behavior of the NICER pipeline is to calculate the pointing offset relative to the source coordinates.</footnote>
<table><caption>**Table 1**
Timing Solution of MAXI J1816</caption><thead><tr><th>Parameter</th><th>Value</th><th>Uncertainty</th></tr></thead><tbody><tr><td>Epoch (MJD)</td><td>59,750</td><td></td></tr><tr><td>P<sub>orb</sub> (d)</td><td>0.20141878</td><td>5 × 10<sup>-8</sup></td></tr><tr><td>a<sub>x</sub> sin i (lt-s)</td><td>0.262949</td><td>1.4 × 10<sup>-5</sup></td></tr><tr><td>T<sub>asc</sub> (MJD)</td><td>59738.8756284</td><td>2.7 × 10<sup>-6</sup></td></tr><tr><td>Eccentricity</td><td>&lt;2 × 10<sup>-4</sup></td><td></td></tr><tr><td colspan="3"><strong>Constant Frequency Model</strong></td></tr><tr><td>ν<sub>0</sub> (Hz)</td><td>528.611105819</td><td>1.0 × 10<sup>-8</sup></td></tr><tr><td>χ<sup>2</sup> / dof</td><td>2167/78</td><td></td></tr><tr><td colspan="3"><strong>Flux-bias Model</strong></td></tr><tr><td>ν<sub>0</sub> (Hz)</td><td>528.611105950</td><td>1.4 × 10<sup>-8</sup></td></tr><tr><td>b (cycles/rate<sup>Γ</sup>)</td><td>4.7 × 10<sup>-5</sup></td><td>5 × 10<sup>-6</sup></td></tr><tr><td>Γ</td><td>1.2</td><td>0.2</td></tr><tr><td>χ<sup>2</sup> / dof</td><td>213/72</td><td></td></tr></tbody></table>

**Note.** The orbital parameters listed are obtained by fitting the constant frequency model to the $t = 5–25$ days data. The orbital parameters obtained with the flux-bias model are consistent with listed values within quoted errors. All reported MJD are barycentric (TDB). Uncertainties are $1\sigma$ errors and upper limits are quoted at 95% confidence.

uncertainty, $A/\sigma_A > 3$. We measured a significant amplitude for the fundamental pulsation in 83 out of 92 pointings, with the nondetections all confined to the final days of the outburst, where the count-rate dropped below the background level. The second harmonic was only required in two pointings near the peak of the outburst. We proceeded by selecting the measured phases of the fundamental and converted them to pulse arrival times. These arrival times were then modeled with TEMPO2 (Hobbs et al. 2006), using a constant pulse frequency and a circular orbital model. We repeated this process of folding the data and fitting the pulse arrival times until the timing solution converged.

The parameters of the obtained timing solution are listed in Table 1, while the resulting pulse amplitudes and phase residuals are shown in Figure 1 (second and third panels). From the figure it is clear that the timing solution describes the decay of the outburst well, but leaves systematic residuals before $t = 5$ days and after $t = 25$ days.

Including a frequency derivative in the model did not meaningfully improve the quality of the fit. Instead, we found that we needed to include terms up to the fifth frequency derivative before the structural phase residuals were reduced. However, such a high-degree polynomial frequency model is plainly unphysical.

In an alternative approach, we attempted to fit the data using the flux-bias model of Bult et al. (2020). We expressed the phase model as

$$ \varphi(t, F) = \varphi_0 + \nu_0 t + \varphi_{\text{orb}}(t) + bF^{\Gamma}, \quad (1) $$

where $\varphi_0$ is an arbitrary phase offset, $\nu_0$ the constant pulsar spin frequency, and $\varphi_{\text{orb}}(t)$ the phase correction associated with the binary orbit. The final term on the right-hand side adds a flux dependent phase bias with scale factor $b$, and power-law index $\Gamma$. We adopted the 0.5–10 keV count-rate as a proportional substitute for flux and set the power-law index to $-1/5$ to model the effect of a phase drift imposed by the moving magnetospheric boundary (Bult et al. 2020). This approach

again only marginally improved the fit. Leaving the power-law index free to vary gives the model sufficient flexibility to account for phase residuals during the outburst rise, but leaves a discrete 0.3 cycle jump in the phase late in the outburst. We list the parameters of this model in Table 1 and show the phase residuals in Figure 1 (fourth panel).

Finally, we investigated the energy dependence of the pulse waveform. We divided the 0.5–10 keV energy range into 50 bins, such that each bin contains a roughly equal number of photons. For each energy bin we folded the entire data set to a pulse profile using the flux-bias timing model and measured the amplitude and phase of the fundamental pulsation. We found that the pulse phase remained constant across all energy bins within measurement uncertainties, while the fractional amplitude increases linearly with energy from $(0.5 \pm 0.2)\%$ at 0.5 keV to $(6.4 \pm 0.1)\%$ at 10 keV.

## 3.2. X-Ray Bursts

We observed 15 type I X-ray bursts from MAXI J1816, all with very similar profiles. To characterize these profiles we manually determined the burst onset, $t_0$, measured the rise time as the time from onset to the peak of the $1/8$ s time resolution light curve, and measured the exponential decay timescale between $[t_0 + 10 \text{ s}, t_0 + 100 \text{ s}]$ (see Table 2). The light curves of bursts #1 through #8 are almost identical, with a rise time of about 8 s and an exponential decay timescale of about 23 s. As the outburst progressed, the X-ray burst profiles showed a modest shift in shape toward faster rise times and shorter decay times. This effect is illustrated in Figure 2.

We investigated the energetics of the X-ray bursts through a time-resolved spectroscopic analysis. First we extracted a preburst spectrum from $[t_0 - 125 \text{ s}, t_0 - 25 \text{ s}]$ (except for burst #4, where we used the final 100 s of the pointing). We then extracted time-resolved spectra from the X-ray bursts by dynamically dividing the burst epochs into multiples of 0.1 s, such that each temporal bin contained at least 1500 events.

We modeled the X-ray burst spectra in the 0.8–10 keV energy range using XSPEC version 12.12 (Arnaud 1996). Subtracting the preburst emission as background, we described each burst spectrum as an absorbed blackbody ($tbabs * bbodyrad$). This model provided statistically acceptable $\chi^2$ values throughout each X-ray burst. The average absorption column density is $(2.36 \pm 0.06) \times 10^{22} \text{ cm}^{-2}$, while the typical peak blackbody temperature and normalization are $1.90 \pm 0.03$ keV and $300 \pm 10 \text{ km}^2$ at 10 kpc, respectively. None of the observed X-ray bursts showed evidence for photospheric radius expansion or spectral lines. In Table 2 we list the bolometric blackbody fluence and peak flux measured in each burst.

While a simple blackbody model is sufficient for individual bursts, we note that the best-fit $\chi^2$ values show a systematic increase around the peak intensity of each X-ray burst. This suggests that the persistent emission may be modestly affected by the burst flux. Indeed, if we fit the spectra at the peak of each X-ray burst jointly with the blackbody temperature and normalization tied across all bursts, then we find that an absorbed blackbody spectrum no longer provides a sufficient description of the data. As an alternative model, we adopt the so-called $f_a$ method (Worpel et al. 2013). We generated a background spectrum for each preburst spectrum using the NICER 3C50 background model (Remillard et al. 2022) and modeled each of the preburst spectra using an absorbed disk
<table><caption>**Table 2**

Detected X-Ray bursts</caption><thead><tr><th>ID</th><th>ObsID</th><th>MJD<br>(TDB)</th><th>Note</th><th>Fluence<br>(×10<sup>-7</sup> erg cm<sup>-2</sup>)</th><th>Peak Flux<br>(×10<sup>-8</sup> erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>Rise<br>(s)</th><th>ε<br>(s)</th><th>τ<br>(s)</th><th>α</th></tr></thead><tbody><tr><td>1</td><td>5533010101</td><td>59738.793937</td><td></td><td>9.5 ± 0.7</td><td>4.1 ± 0.3</td><td>7.6</td><td>23.2 ± 0.2</td><td>23.1</td><td></td></tr><tr><td>2</td><td>5533010101</td><td>59738.976061</td><td>SAA</td><td>9.9 ± 0.8</td><td>4.6 ± 0.4</td><td>8.8</td><td>24.3 ± 0.2</td><td>21.7</td><td></td></tr><tr><td>3</td><td>5533010102</td><td>59739.366594</td><td></td><td>9.7 ± 0.7</td><td>4.2 ± 0.4</td><td>7.4</td><td>23.5 ± 0.3</td><td>23.1</td><td></td></tr><tr><td>4</td><td>5533010103</td><td>59740.537525</td><td>Tail</td><td>>6.7 ± 0.5</td><td>>4.0 ± 0.4</td><td>...</td><td>23.6 ± 0.3</td><td>...</td><td></td></tr><tr><td>5</td><td>5533010103</td><td>59740.790575</td><td></td><td>9.1 ± 0.7</td><td>4.3 ± 0.4</td><td>6.9</td><td>23.5 ± 0.3</td><td>21.2</td><td></td></tr><tr><td>6</td><td>5533010104</td><td>59741.633451</td><td></td><td>9.4 ± 0.7</td><td>4.3 ± 0.3</td><td>9.4</td><td>22.7 ± 0.2</td><td>21.9</td><td></td></tr><tr><td>7</td><td>5533010105</td><td>59742.413764</td><td>SAA</td><td>9.5 ± 0.7</td><td>4.3 ± 0.3</td><td>7.9</td><td>23.6 ± 0.3</td><td>22.3</td><td></td></tr><tr><td>8</td><td>5533010105</td><td>59742.470702</td><td></td><td>9.4 ± 0.7</td><td>4.2 ± 0.3</td><td>7.3</td><td>22.8 ± 0.2</td><td>22.4</td><td>45 ± 3</td></tr><tr><td>9</td><td>5533010106</td><td>59743.103625</td><td>Rise</td><td>>0.8 ± 0.1</td><td>>3.1 ± 0.2</td><td>>4.7</td><td>...</td><td>...</td><td></td></tr><tr><td>10</td><td>5533010801</td><td>59746.344881</td><td>SAA</td><td>9.3 ± 0.6</td><td>4.6 ± 0.4</td><td>5.6</td><td>21.4 ± 0.2</td><td>20.3</td><td></td></tr><tr><td>11</td><td>5533010901</td><td>59747.308052</td><td></td><td>9.0 ± 0.6</td><td>4.1 ± 0.3</td><td>5.7</td><td>20.6 ± 0.2</td><td>21.7</td><td></td></tr><tr><td>12</td><td>5533010901</td><td>59747.567004</td><td></td><td>8.8 ± 0.6</td><td>3.9 ± 0.3</td><td>6.6</td><td>19.4 ± 0.3</td><td>22.3</td><td></td></tr><tr><td>13</td><td>5533011301</td><td>59751.757114</td><td></td><td>8.4 ± 0.6</td><td>4.4 ± 0.3</td><td>4.8</td><td>19.5 ± 0.2</td><td>19.1</td><td></td></tr><tr><td>14</td><td>5533011601</td><td>59754.160772</td><td></td><td>8.3 ± 0.5</td><td>4.2 ± 0.3</td><td>6.6</td><td>18.6 ± 0.2</td><td>19.5</td><td>58 ± 4</td></tr><tr><td>15</td><td>5533011601</td><td>59754.288382</td><td></td><td>8.1 ± 0.5</td><td>4.0 ± 0.3</td><td>4.8</td><td>19.3 ± 0.2</td><td>20.2</td><td>61 ± 5</td></tr></tbody></table>

Note. The MJD column lists the burst onset times. Fluxes are bolometric. Columns ε and τ give the exponential decay timescale and the ratio burst fluence to peak flux, respectively. The α measurement of burst 14 follows from the 3 hr recurrence time obtained from simultaneous NuSTAR coverage (see the text). All uncertainties give 1σ errors.

blackbody plus a thermally Comptonized continuum (nthcomp; Zdziarski et al. 1996; Życki et al. 1999). Finally, the X-ray burst spectra were modeled using an absorbed blackbody plus the preburst model, where the preburst component was multiplied with a variable factor, $f_a$. This approach improved the $\chi^2$/ dof from 985/863 to 945/862 (factor 10 improvement in the p-value) with $f_a = 1.22 \pm 0.04$ at peak burst intensity.

### 3.2.1. Burst Recurrence Time

Dividing the total unfiltered exposure (111 ks) by the number of detected X-ray bursts, we estimate the average burst recurrence time at 2.1 hr. This average is close to the actual spacing observed between X-ray bursts: we find bursts #7 and #8 are separated by 1.4 hr, while bursts #14 and #15 are separated by 3.1 hr.

MAXIJ1816 was also observed with NuSTAR on 2022 June 23 (MJD 59753; Chauhan et al. 2022). This NuSTAR observation contains four X-ray bursts, the first two of which were not observed with NICER and the latter two being #14 and #15 in our sample. This train of four X-ray bursts observed with NuSTAR is consistent with a regular recurrence time of 3.0 hr, while a shorter recurrence time for these bursts (1.5 hr) is ruled out by the joint NICER and NuSTAR coverage. Hence, these results suggest that MAXIJ1816 is a regular burster with a recurrence time that lengthens over the course of the outburst.

The α factor is defined as the ratio of the persistent fluence over the burst fluence (Galloway et al. 2020), and can be estimated as

$$ \alpha = \frac{F_{\text{persist}} \Delta t}{E_{\text{burst}}} \qquad (2) $$

where $F_{\text{persist}}$ is the bolometric flux of the persistent emission, $\Delta t$ the time between successive X-ray bursts, and $E_{\text{burst}}$ the fluence of the burst. We estimate the $F_{\text{persist}}$ by adding a cflux component to the persistent spectra and measuring the flux between 0.01 and 100 keV, while we obtain the burst

fluence from the time-resolved spectroscopy. For the three X-ray bursts with a reliable measurement of the recurrence time (#8, #14, and #15) we find α values of 46, 58, and 61 (respectively).

### 3.2.2. Burst Oscillations

To search the X-ray bursts for the presence of coherent burst oscillations, we employed a sliding window search method. For each X-ray burst, we constructed a light curve at a 1/8192 s time resolution and applied a window to this light curve of the duration $T=2, 4, 8$ s. We then moved the window across the light curve in steps of $T/4$. At each window position we calculated the power density spectrum and searched for single bin powers between $528.6 \pm 5$ Hz that exceeded the $3\sigma$ detection threshold calculated from the counting noise distribution (correcting for the number of trials; the number of windows times the number of powers per window). We applied this search strategy using all events in the 0.5–10 keV energy range. No burst oscillation candidates were found.

## 4. Discussion

We have presented the discovery of 528.6 Hz pulsations from MAXIJ1816. Through a timing analysis of the pulsations, we measured the binary ephemeris reported in Table 1. From this ephemeris we find that the pulsar mass function is $f_x = 4.8 \times 10^{-3} M_{\odot}$, implying a minimum companion mass of $M_2 > 0.10 M_{\odot}$ for a canonical $1.4 M_{\odot}$ neutron star.

If we assume that the companion star fills its Roche lobe, we can use the Roche lobe radius (Eggleton 1983) to calculate both the mass and radius of the companion star as a function of the binary inclination. This empirical mass–radius relation intersects with the theoretical mass–radius relation for a zero-age main-sequence star (Tout et al. 1996) at a companion mass of $0.55 M_{\odot}$ (Figure 3). If the companion star is evolved or ablated by the irradiation of the compact object, however, it will tend toward a larger radius for the same stellar mass, which means that we can treat this intersection as an upper limit on the
<figure><img src="image_4.png" /><figcaption>Figure 2. Light curves of the 15 X-ray bursts observed from MAXI J1816. These light curves are calculated in the 0.5–10 keV energy range at 1 s resolution, with the preburst count-rate subtracted. In each panel a different subset of the bursts are highlighted in color (with the burst number as in Table 2), while the remaining bursts are shown in gray for comparison.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Mass–radius relation of the companion star from the binary ephemeris (solid line) compared to the theoretical mass–radius relation of a zero-age main-sequence star (dotted line).</figcaption></figure>

companion mass. Indirectly, this mass limit then implies a lower limit on the binary inclination of 13°.

We found that the pulse phase shows a complex evolution with time. The greater part of the observed outburst, between t = 5 and 25 d, could be modeled using a constant pulse frequency model. Outside this time range the pulse phase shows systematic residuals relative to the model. Similar structural patterns in the phase residuals have been observed in a number of other AMXPs (Burderi et al. 2007; Hartman et al. 2008; Patruno et al. 2009; Bult et al. 2020; Sanna et al. 2020). The common explanation for this behavior is that either the accretion torque is measurably changing the stellar spin frequency over the course of the outburst, or that the hotspot position on the stellar surface is not fixed, but shifts in response to the changing accretion geometry. Of course these effects are not exclusive, and could both be contributing factors.

For MAXI J1816 we were able to account for the phase residuals observed during the outburst rise by employing a flux-bias model for the pulsar phase. Interestingly, this model prefers a power-law index of Γ = 1.2, which is a shallower dependence than the Γ ≳ 2 expected from a phase drift imposed through a cumulative accretion torque (Sanna et al. 2020). Yet, this index is not consistent with the hotspot position changes

predicted from numerical simulations of accreting pulsars either (Kulkarni & Romanova 2013). Hence, if the hotspot is moving, its phase bias is not driven by the changing radius of the magnetosphere, but is instead nearly proportional to intensity and thus the mass accretion rate.

Even with the flux-bias model, we could not account for a discrete 0.3 cycle jump in the pulse phase observed in the final day of the outburst. Similar late-time phase jumps have been observed from SAX J1808.4–3658 (Burderi et al. 2006; Hartman et al. 2008), and have been attributed to a viewing geometry driven by the receding accretion disk (Ibragimov & Poutanen 2009). Given that the phase jump in MAXI J1816 occurs right before the source drops below the detection level, it is plausible that a similar mechanism may be at play here.

## 4.1. Constraints from the X-Ray Bursts

We found that all observed X-ray bursts showed a very similar evolution, both in the light curve and in their spectroscopy. The burst durations and α values point to ignition in a hydrogen-rich environment (Lewin et al. 1993; Galloway & Keek 2021), indicating that the accreted material and thus the donor star must be hydrogen-rich. The peak fluxes are consistent within errors, with an average of (4.3 ± 0.1) × 10⁻⁸ erg s⁻¹ cm⁻². Depending on the hydrogen abundance in the burst fuel, the expected Eddington luminosity is 2.2–3.8 × 10³⁸ erg s⁻¹ (Kuulkers et al. 2003), which yields an upper limit on the distance of d < 6.5–8.6 kpc. The lower end of this range is associated with hydrogen-rich bursts, and is therefore preferred.

For three of the observed X-ray bursts we obtained a reliable measurement for α, finding that this factor increases from 46 for burst #8 to about 60 for bursts #14 and #15. Such evolution usually indicates a shift in the average hydrogen abundance of the burst fuel (Galloway et al. 2004), which would be consistent with the smaller fluence and decay timescale seen in later bursts. Specifically, if we assume that all accreted matter burns during an X-ray burst, then α follows from theory as (Galloway et al. 2020)

$$ \alpha = \frac{Q_{\text{gravity}}}{Q_{\text{nuclear}}} (1 + z) \frac{\xi_{\text{burst}}}{\xi_{\text{disk}}} \quad (3) $$
where $Q_{\text{gravity}} = GM_{\text{ns}}/R_{\text{ns}}$ is the gravitational binding energy, $Q_{\text{nuclear}} = 1.35 + 6.05\bar{X}$ MeV nucleon$^{-1}$ the nuclear energy released by accreted matter (Goodwin et al. 2019; with $\bar{X}$ the average hydrogen abundance in the burst fuel), $z$ the gravitational redshift, and $\xi_{\text{burst}}/\xi_{\text{disk}}$ the ratio of the burst and disk anisotropy (Fujimoto 1988). If the CNO cycle is stably burning hydrogen between the X-ray bursts then shorter recurrence times give a higher $\bar{X}$ and thus a higher energy release per accreted nucleon (Galloway & Keek 2021). For solar abundances the observed increase in recurrence time would then increase $\alpha$ by about 10%, well short of what is needed to explain the observations. This suggests that some additional physical process is gradually changing the ignition condition either as a function of time or mass accretion rate. Whether that missing piece relates to an evolving accretion geometry, an inertia in the burst train (Woosley et al. 2004; Johnston et al. 2018), a temperature evolution in the neutron star crust (Chenevez et al. 2011), or perhaps some mixing effect related to the ignition latitude (Cavecchi et al. 2020) remains an open question.

This work made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC). P.B. acknowledges support from NASA through the NICER Guest Observer Program and the CRESST II cooperative agreement (80GSFC21M0002). S.G. acknowledges the support of the CNES. D.A. acknowledges support from the Royal Society. G.C.M. was partially supported by Proyecto de Investigación Plurianual (PIP) 0102 (Nacional de Investigaciones Científicas y Técnicas (CONICET)) and from PICT-2017-2865 (Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT)). NICER science team members at NRL are supported by NASA under Interagency Agreement NNG200808A.

Facilities: ADS, HEASARC, NICER.

Software: heasoft (v6.30), nicerdas (v9), xspec (Arnaud 1996), tempo2 (Hobbs et al. 2006).

## ORCID iDs

Peter Bult https://orcid.org/0000-0002-7252-0991

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

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

Jérôme Chenevez https://orcid.org/0000-0002-4397-8370

Elizabeth C. Ferrara https://orcid.org/0000-0001-7828-7708

Keith C. Gendreau https://orcid.org/0000-0001-7115-2819

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

Tolga Güver https://orcid.org/0000-0002-3531-9842

Wataru Iwakiri https://orcid.org/0000-0002-0207-9010

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

Giulio C. Mancuso https://orcid.org/0000-0001-9822-6937

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

Mason Ng https://orcid.org/0000-0002-0940-6563

Andrea Sanna https://orcid.org/0000-0002-0118-2649

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

Zorawar Wadiasingh https://orcid.org/0000-0002-9249-0515

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

## References

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

Beauchamp, I., Belvedere, S., Hernandez, M., et al. 2022, ATel, 15481, 1

Bildsten, L. 1998, ApJL, 501, L89

Bright, J., Russell, T., Evangelia, T., et al. 2022, ATel, 15484, 1

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

Bult, P., Chakrabarty, D., Arzoumanian, Z., et al. 2020, ApJ, 898, 38

Bult, P. M., Ng, M., Altamirano, W. I. D., et al. 2022a, ATel, 15425, 1

Bult, P. M., Sanna, A., Ng, M., et al. 2022b, ATel, 15431, 1

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

Burderi, L., Di Salvo, T., Lavagetto, G., et al. 2007, ApJ, 657, 961

Cavecchi, Y., Galloway, D. K., Goodwin, A. J., Johnston, Z., & Heger, A. 2020, MNRAS, 499, 2148

Chauhan, J., Lohfink, A., Bharali, P., et al. 2022, ATel, 15470, 1

Chenevez, J., Altamirano, D., Galloway, D. K., et al. 2011, MNRAS, 410, 179 de Martino, D., D'Avanzo, P., Ambrosino, F., et al. 2022, ATel, 15479, 1

Di Salvo, T., & Sanna, A. 2022, in Millisecond Pulsars, ed. S. Bhattacharyya, A. Papitto, & D. Bhattacharya (Cham: Springer), 87

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

Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, IPNPR, 42, 1

Fujimoto, M. Y. 1988, ApJ, 324, 995

Galloway, D. K., Cumming, A., Kuulkers, E., et al. 2004, ApJ, 601, 466

Galloway, D. K., & Keek, L. 2021, in Timing Neutron Stars: Pulsations, Oscillations and Explosions, ed. T. M. Belloni, M. Méndez, & C. Zhang (Belin: Springer), 209

Galloway, D. K., in't Zand, J., Chenevez, J., et al. 2020, ApJS, 249, 32

Goodwin, A. J., Heger, A., & Galloway, D. K. 2019, ApJ, 870, 64

Hameury, J. M. 2020, AdSpR, 66, 1004

Hartman, J. M., Patruno, A., Chakrabarty, D., et al. 2008, ApJ, 675, 1468

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

Ibragimov, A., & Poutanen, J. 2009, MNRAS, 400, 492

Johnston, Z., Heger, A., & Galloway, D. K. 2018, MNRAS, 477, 2112

Kennea, J. A., Evans, P. A., & Negoro, H. 2022a, ATel, 15421, 1

Kennea, J. A., Evans, P. A., & Negoro, H. 2022b, ATel, 15467, 1

Kulkarni, A. K., & Romanova, M. M. 2013, MNRAS, 433, 3048

Kuulkers, E., den Hartog, P. R., in't Zand, J. J. M., et al. 2003, A&A, 399, 663

Lasota, J.-P. 2001, NewAR, 45, 449

Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1993, SSRv, 62, 223

Negoro, H., Serino, M., Iwakiri, W., et al. 2022, ATel, 15418, 1

Patruno, A., & Watts, A. L. 2021, in Timing Neutron Stars: Pulsations, Oscillations and Explosions, ed. T. M. Belloni, M. Méndez, & C. Zhang (Berlin: Springer), 143

Patruno, A., Wijnands, R., & van der Klis, M. 2009, ApJL, 698, L60

Poutanen, J., & Gierlinski, M. 2003, MNRAS, 343, 1301

Psalis, D., & Chakrabarty, D. 1999, ApJ, 521, 332

Remillard, R. A., Loewenstein, M., Steiner, J. F., et al. 2022, AJ, 163, 130

Sanna, A., Burderi, L., Gendreau, K. C., et al. 2020, MNRAS, 495, 1641

Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257

Wijnands, R., & van der Klis, M. 1998, ApJL, 507, L63

Woosley, S. E., Heger, A., Cumming, A., et al. 2004, ApJS, 151, 75

Worpel, H., Galloway, D. K., & Price, D. J. 2013, ApJ, 772, 94

Zdziarski, A. A., Johnson, W. N., & Magdziarz, P. 1996, MNRAS, 283, 193

Życki, P. T., Done, C., & Smith, D. A. 1999, MNRAS, 309, 561