# Long-term Study of the 2020 Magnetar-like Outburst of the Young Pulsar PSR J1846-0258 in Kes 75

R. Sathyaprakash<sup>1,2,3</sup>, N. Rea<sup>1,2</sup>, F. Coti Zelati<sup>1,2</sup>, A. Borghese<sup>4,5</sup>, M. Pilia<sup>6</sup>, M. Trudu<sup>6,7</sup>, M. Burgay<sup>6</sup>, R. Turolla<sup>8,9</sup>, S. Zane<sup>9</sup>, P. Esposito<sup>3,10</sup>, S. Mereghetti<sup>10</sup>, S. Campana<sup>11</sup>, D. Götz<sup>12</sup>, A. Y. Ibrahim<sup>1,2</sup>, G. L. Israel<sup>13</sup>, A. Possenti<sup>6</sup>, and A. Tiengo<sup>3,10</sup>

<sup>1</sup> Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans s/n, 08193, Barcelona, Spain; rajath.sathyaprakash@iusspavia.it

<sup>2</sup> Institut d'Estudis Espacials de Catalunya (IEEC), Esteve Terradas 1, 08860, Castelldefels, Barcelona, Spain

<sup>3</sup> Scuola Universitaria Superiore IUSS Pavia, Palazzo del Broletto, piazza della Vittoria 15, I-27100 Pavia, Italy

<sup>4</sup> Instituto de Astrofísica de Canarias (IAC), Vía Láctea s/n, La Laguna 38205, S/C de Tenerife, Spain

<sup>5</sup> European Space Astronomy Centre (ESA/ESAC), Villanueva de la Cañada (Madrid), Spain

<sup>6</sup> INAF Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius, Italy

<sup>7</sup> Università degli Studi di Cagliari, Dipartimento di Fisica, SP Monserrato-Sestu km 0.7, I-09042 Monserrato (CA), Italy

<sup>8</sup> Dipartimento di Fisica e Astronomia “Galileo Galilei,” Università di Padova, via F. Marzolo 8, I-35131 Padova, Italy

<sup>9</sup> Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey, RH5 6NT, UK

<sup>10</sup> INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica di Milano, via A. Corti 12, I-20133 Milano, Italy

<sup>11</sup> INAF-Osservatorio Astronomico di Brera, Via Bianchi 46, Merate (LC), I-23807, Italy

<sup>12</sup> AIM-CEA/DRF/Irfu/Département d'Astrophysique, CNRS, Université Paris-Saclay, Université de Paris Cité, Orme des Merisiers, F-91191 Gif-sur-Yvette, France

<sup>13</sup> INAF-Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone, Italy
Received 2024 January 12; revised 2024 September 27; accepted 2024 September 27; published 2024 November 13

## Abstract

Magnetar-like activity has been observed in a large variety of neutron stars. PSR J1846−0258 is a young 327 ms radio-quiet pulsar with a large rotational power (~8 × 10<sup>36</sup> erg s<sup>−1</sup>), and resides at the center of the supernova remnant Kes 75. It is one of the rare examples of a high-magnetic-field pulsar showing characteristics both of magnetars and radio pulsars, and can thus provide important clues on the differences in the emission mechanisms between these two classes. In 2006, PSR J1846−0258 was detected to undergo an outburst for the first time, accompanied by a large flux increase, millisecond X-ray bursts, significant spectral changes, and a large timing glitch. In the period between 2020 May and June, after 14 yr of quiescent stable emission, the source underwent a second magnetar-like outburst, which was followed up with several observations by Neutron Star Interior Composition Explorer, XMM-Newton, NuSTAR, and Swift. In this work, we report on the long-term timing and X-ray spectral properties of the source following the 2020 outburst, and place upper limits on any source activity at radio wavelengths. We demonstrate that the pulsed flux increased by a factor >6 during the outburst, followed by nontrivial variability in the spin-down rate. Our timing analysis shows that the spin frequency and its derivative are clearly affected by magnetospheric activity due to the outburst. We find hints for an oscillation in the frequency derivative with a timescale of 50–60 days, recovering later on to stable quiescence.

*Unified Astronomy Thesaurus concepts: Magnetars (992); X-ray transient sources (1852); Rotation powered pulsars (1408); Neutron stars (1108)*

*Materials only available in the online version of record: data behind figure, machine-readable table*

## 1. Introduction

Isolated neutron stars (NSs) are the compact remnants of core-collapse supernova explosions, characterized by fast rotation velocities and large magnetic fields. Observational surveys limited to the Milky Way alone have discovered more than three thousand isolated NSs in the past few decades. A handful of around 30 such objects, known as magnetars, were historically observed as flaring X-ray and gamma-ray sources and are believed to be powered by extreme magnetic fields (i.e., ~10<sup>14</sup>–10<sup>15</sup> G; for a detailed review, see, e.g., P. Esposito et al. 2021 and V. M. Kaspi & A. M. Beloborodov 2017) with relatively slow spin periods of 1–10 s and spin-down ages ranging from 10<sup>2</sup> to 10<sup>4</sup> yr. Among the key attributes of magnetars is that their X-ray luminosity appears to far exceed

the available spin-down power, especially during periods of short-duration (i.e., millisecond) bursts and long-duration outbursts. Their spectra during quiescence are typically modeled by blackbody emission with temperatures ~0.3–1.0 keV, with a hard X-ray tail.

The difference between the emission of rotation-powered pulsars (RPPs) and magnetars has been studied for many years. Two decades ago, the former were observed mainly as stable radio emitters, while the latter were observed as radio-quiet bright X-ray sources. This was attributed to a distinct emission reservoir and physics between the two classes. The discovery of radio-emitting magnetars (F. Camilo et al. 2006), or magnetar-like emission in allegedly canonical radio pulsars (e.g., F. P. Gavriil et al. 2008; R. F. Archibald et al. 2016) and low magnetic field magnetars (e.g., N. Rea et al. 2010), has completely changed our understanding of the magnetic-powered emission across the isolated NS zoo. Indeed, more recent studies show that the population of isolated magnetars appears to not only be restricted to Galactic sources, and those

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
that have been observed to emit millisecond duration radio flares have been posited to be excellent candidates for fast radio bursts (FRBs; e.g., F. Y. Wang et al. 2022), a few of which have been deduced to be extragalactic in origin due to their large dispersion measures.

Located within the core of the supernova remnant (SNR) Kes 75, PSR J1846−0258 (first discovered by E. V. Gotthelf et al. 2000) is one of the youngest pulsars in our Galaxy with a characteristic age ~700 yr, a spin period of 327 ms, and a distance of 5.8 kpc (J. P. W. Verbiest et al. 2012).<fnref n="14" /> It has a dipolar magnetic field of $B \sim 4.5 \times 10^{13}$ G and has historically been radio-quiet but a very powerful X-ray emitter. Similar to other RPPs, it has a large rotational power ($\dot{E} \sim 8 \times 10^{36}$ erg s⁻¹) and lies at the center of a bright pulsar wind nebula (PWN).

For many years prior to 2006, it behaved as a very stable X-ray pulsar. However, in 2006, the source displayed a distinctly magnetar-like behavior, undergoing a dramatic increase in pulsed flux that coincided with the detection of X-ray bursts lasting a few milliseconds (F. P. Gavriil et al. 2008), and showed significant spectral changes (H. S. Kumar & S. Safi-Harb 2008) and a large timing glitch (L. Kuiper & W. Hermsen 2009a). Therefore, this object has been recognized as an interesting hybrid between magnetars and rotation-powered pulsars (RPPs) that helps to understand the links between these two different classes of NSs.

In 2009, a Chandra observation caught the source in quiescence (M. A. Livingstone et al. 2011), with a braking index ($n = 2.16 \pm 0.13$) that differed significantly from that measured during the pre-outburst epoch ($n = 2.65 \pm 0.01$) in 2006 (M. A. Livingstone et al. 2006). Indeed, further investigations showed that the braking index had changed permanently following the magnetar-like outburst and its radiatively loud glitch (R. F. Archibald et al. 2015). Thus far, no physical model can consistently explain the long-term decrease in the braking index.

Other works have focused on studying the emission of the pulsar and its wind nebula using XMM-Newton, NuSTAR, and Chandra observations (e.g., E. V. Gotthelf et al. 2021). By performing phase-resolved spectroscopy, E. V. Gotthelf et al. (2021) found that the joint XMM-Newton and NuSTAR spectrum of the pulsar is characterized by a power-law model in the 2–50 keV energy band. However, when analyzing a broader energy range (2 keV–100 MeV) using data from RXTE, INTEGRAL, and Fermi, L. Kuiper et al. (2018) found evidence for the presence of spectral curvature in the pulsed emission of PSR J1846−0258. This phenomenology could potentially be explained in the context of the Outer Gap, Slot Gap, or Polar Cap models, in which nonthermal gamma-ray photons are emitted due to the synchrotron and inverse Compton processes by accelerating electrons. L. Kuiper et al. (2018) also reported weak pulsed emission in the soft gamma-ray range (30–100 MeV) analyzing Fermi-LAT data collected during an ~8 yr period, after the source had returned to a quiescent state following its 2006 outburst. In more recent work, S. M. Straal et al. (2023) used Fermi-LAT data to observe the spectral energy distribution of the source beyond 100 MeV, and thereby placed constraints on the properties of the PWN (i.e., its magnetization and particle energy spectrum),

in addition to the birth properties of the NS and supernova ejecta.

In the period between 2020 May and June, the source underwent a second outburst (followed by a short 100 ms duration burst in 2020 August; see H. A. Krimm et al. 2020 and M. Uzuner et al. 2023) after a long quiescent period that had lasted for more than 14 yr (H. A. Krimm et al. 2020; S. Laha et al. 2020). This event triggered the onset of several observations with the Neutron Star Interior Composition Explorer (NICER) and the Neil Gehrels Swift Observatory (Swift). These observations revealed the occurrence of a spin-up glitch (L. Kuiper et al. 2020) accompanied by a significant increase in pulsed flux, reaching over 5 times the quiescent level, with the pulsed spectrum characterized by the combination of a blackbody and a power-law component between 1 and 70 keV (C.-P. Hu et al. 2023). Spectral analysis indicated that the flux increase was primarily caused by the emergence of the thermal component, which gradually faded away, while the power-law component remained constant throughout the outburst (C.-P. Hu et al. 2023). Observations with the Deep Space Network (DSN) 34 m radio instrument early in the outburst resulted in a nondetection of radio pulsations from the source, consistent with previous observations (see W. A. Majid et al. 2020; H. Blumer et al. 2021 and M. B. Mickaliger et al. 2020).

In the present work, we report on the long-term timing and spectral properties of the source following the recent outburst using four NuSTAR observations, three XMM-Newton observations, and eight Swift observations, all triggered in the aftermath of the millisecond burst (see H. A. Krimm et al. 2020). We also include several NICER observations obtained over 3 yr from 2018 August until 2021 August (see also C.-P. Hu et al. 2023).

## 2. Data Reduction

The X-ray observations analyzed in this work are listed at the end of this article. The EPIC-pn XMM-Newton data were acquired in large window mode, while the EPIC-MOS data were collected in full-frame mode. Owing to their insufficient timing resolution (2.6 s), data from the MOS detectors were excluded for timing analysis. To process the XMM-Newton data, we used the Science Analysis Software (SAS v.18). We extracted a light curve between 10 and 12 keV for each observation (in the entire field of view) using `evselect` to identify and filter out periods of intense background flares. We then generated spatially filtered EPIC-pn event files, extracting source events within a circular region of 30" in the energy range of 1–12 keV, centered on the most precise available Chandra source position (R.A. = 281.6039167, decl. = −2.9750278; D. J. Helfand et al. 2003). The same position was adopted for all the other instruments. Background events were accumulated in a neighboring circular region of 30" radius.

The NuSTAR data were processed with the `nupipeline` script and source events were collected from a circular region of 50" radius, within the recommended energy range of 3–79 keV. Background photons were extracted from a neighboring source-free circular region of the same size. The `nupipeline` script excludes passages of the spacecraft through the South Atlantic Anomaly (SAA), with a recently updated algorithm that simultaneously monitors the count rates from the high gain shield and a more efficient cadmium zinc telluride detector, providing a cleaned event file that is suitable for timing analysis.

<footnote n="14">¹⁴ Distance derived from analyzing H I and ¹³CO maps of the Kes 75 remnant (D. A. Leahy & W. W. Tian 2008).</footnote>
The NICER data were calibrated with the nicer12 pipeline using the heasoft package (version 6.31; NASA High Energy Astrophysics Science Archive Research Center (HEA-SARC) 2014), which excludes time intervals within the SAA, in addition to overshoot and undershoot events caused by the charged particle background and solar optical photons interacting with the detector. Information from each Measurement Power Unit was combined to produce a cleaned event file.

Finally, the Swift-XRT data were acquired in windowed timing mode (WT; timing resolution of 1.8 ms). We reprocessed the data adopting standard cleaning criteria with the task xrtpipeline. The source photons were accumulated from a 10 pixel circular region (1 pixel = 2"/36), while an annulus with inner and outer radii of 60 and 100 pixels, respectively, was used to extract the background events.

We barycentered event files from all instruments using the barycorr tool available as part of the heasoft package (v6.31), with the source position adopted from the Chandra observations (see above) and the latest DE-430 solar system ephemeris. For the XMM-Newton data, we used the barycen tool to perform the barycentric corrections. For the reasons discussed in Section 3.3, we refer to all observations after 2020 July 25 (MJD 59,055) as the post-outburst period, and observations before this date as the pre-outburst period.

### 3. Timing Analysis

During its 2006 outburst, PSR J1846−0258 experienced a spin-up glitch event (L. Kuiper & W. Hermsen 2009b; M. A. Livingstone et al. 2011) that coincided with a pulsed flux increase (F. P. Gavriil et al. 2008) and the detection of four millisecond X-ray bursts. This was followed by a period of overrecovery, such that the glitch could not be described by a simple exponentially decaying model (M. A. Livingstone et al. 2011). The magnitude of the overrecovery was larger than ever observed in other RPPs or magnetars, indicating that the model of vortex unpinning in the superfluid crust of the NS invoked to explain standard glitches in RPPs (see, e.g., B. Haskell & A. Melatos 2015) likely does not hold for PSR J1846−0258 (see M. A. Livingstone et al. 2011 for more details). Here, we describe our timing analysis of PSR J1846−0258 throughout its 2020 outburst to monitor its spin evolution.

#### 3.1. Phase-coherent Timing

To measure the spin parameters of PSR J1846−0258 we attempted to derive a coherent timing solution using the standard phase-connection approach carried out with pulsar timing software. We first assigned a rotational phase to each photon based on a provisional ephemeris using the photo-nphase tool of the PINT software package (J. Luo et al. 2021). The events for each observation were folded to produce integrated pulse profiles, with the adopted energy range discussed in the previous section. The pulse profile that had the highest H-test score was then fitted with a multi-Gaussian model to produce a standard profile template. We then used the photon_toa.py tool of the NICERsoft package to extract as many reliable barycentered pulse time of arrivals (TOAs) as possible within each observation to monitor better the pulsar rotations during the phase-connection process.<fnref n="15" /> Initially, the only fitting parameters of the timing model were the spin

frequency $v$, the first frequency derivative $\dot{v}$, and a phase offset, or “jump,” for every gap between groups of TOAs from the same observation.

We then proceeded with the phase connection with the TEMPO pulsar timing software (D. Nice et al. 2015) by determining the rotation count between groups of TOAs (i.e., eliminating jump parameters) recursively from the shortest to longest gaps. The presence of strong timing noise in PSR J1846−0258 during its outburst required that we include the spin-frequency second derivative $\ddot{v}$ in the model early in the connection process in order to maintain coherence. As we reached the longest gaps, the rotation count became ambiguous: multiple statistically acceptable (reduced $\chi^2 < 2$) and distinct solutions were found when connecting TOAs across these long gaps.

This was problematic in particular during the first 2 months of the outburst (2020 June and July)—only three NICER observations of PSR J1846−0258 were carried out, and we find evidence of (unconstrained) timing anomalies as reported by C.-P. Hu et al. (2023). We attempted to eliminate phase ambiguities using the Dracula algorithm (P. C. C. Freire & A. Ridolfi 2018), but were unable to converge to a single global ephemeris.

To better monitor the spin-down behavior during the outburst, we instead opt to measure local phase-connected timing solutions with TEMPO. TOAs from consecutive observations spanning only a short time range of a few days were used for each such that the solutions were only sensitive to the $v$ and $\dot{v}$ parameters. Local solutions were produced from data taken approximately 5 months prior to the onset of the outburst (i.e., 2020 June–August) until one year after the outburst. Our best-fit $v$ and $\dot{v}$ measurements are shown as black points in Figure 2.

A timing analysis of PSR J1846−0258 during its 2020 outburst was also carried out by C.-P. Hu et al. (2023) using NICER-XTI and Swift-XRT monitoring data—their results are also shown in Figure 2 (solid blue line). We compare our measurements of the spin frequency and spin-down rate with the solutions from C.-P. Hu et al. (2023), which are generally in agreement with their values excluding a few observations (see Figure 2). C.-P. Hu et al. (2023) detected a strong glitch event of fractional size $\Delta v/v \sim 3 \times 10^{-6}$ early during the outburst, between the epochs MJD 58,996 and 59,055. Similarly to C.-P. Hu et al. (2023), we detect a timing irregularity that is consistent with a glitch of estimated magnitude $\Delta v/v \sim 2.5-4 \times 10^{-6}$ during the first month of the outburst (MJD 58,996–59,026).

In summary, we find PSR J1846−0258 experienced rapid torque fluctuations during its outburst. Prior to 2020 June, the spin-down rate was stable at $\dot{v} \sim -6.6 \times 10^{-11}$ Hz s⁻¹. A large $\sim 10^{-5}$ Hz increase in $v$ and a 3% decrease in $\dot{v}$ were measured in the 2020 June 26 (MJD 59,026) NICER data, while the spin measurements the following month (i.e., 2020 July) appeared to have evolved toward values more similar to the pre-outburst $v$ and $\dot{v}$. Such behavior in spin and spin-down rate is consistent with the expected spin evolution of a pulsar experiencing a spin-up glitch event (see, e.g., G. Hobbs et al. 2010) between 2020 May 27 and 2020 June 26. The sparse sampling in this period yields poor constraints on the glitch epoch and thus several statistically valid timing models corresponding to different fractional frequency jumps $\Delta v/v$ in the range of $\sim 2-5 \times 10^{-6}$—of similar size to the glitch event

<footnote n="15">¹⁵ https://github.com/paulray/NICERsoft</footnote>
that accompanied the 2006 outburst of PSR J1846−0258 by L. Kuiper & W. Hermsen (2009a).

Between 2020 August 2 and 21 (MJDs 59,063–59,082), unlike typical post-glitch recovery, the spin-down rate steadily decreased until it reached a minimum of ~−7.1 × 10⁻¹¹ Hz s⁻¹. Over the following weeks $\dot{\nu}$ increased again but at an irregular rate. Our results, shown in Figure 2, are consistent with the spin-down rates reported by C.-P. Hu et al. (2023). Due to visibility constraints, the monitoring of PSR J1846−0258 with NICER was interrupted at the end of 2020 November. When observations resumed in 2021 March, the pulsar’s torque was stable again, but with a spin-down rate $\dot{\nu} \sim -6.7 \times 10^{-11}$ Hz s⁻¹, slightly lower than its pre-outburst value.

## 3.2. Gaussian Process Regression

The timing residuals corresponding to the post-outburst period (between MJD 59,055–59,171), and extracted in the manner described in the previous section display a complex quasi-periodic structure superposed on the secular spin-down of the source. In order to model these residuals, we implemented Gaussian process regression (GPR) via the scikit-learn library (F. Pedregosa et al. 2011). GPR is a nonparametric regression method that is useful in cases where it is not possible to make clear assumptions about a physical model that the data might satisfy. Given a set of data points $(x, y)$, a GPR model can fit (and interpolate between) these data points with functions sampled from a multivariate Gaussian distribution. These functions are characterized by a mean vector and a covariance matrix $k(x, x')$, whose diagonal elements describe the variance in the dependent variable $y$ being modeled, and the off-diagonal elements describe the similarity of the function (or its shape) between any two points $(y(x), y(x'))$. There are several functional forms available for $k(x, x')$, and the choice of which to use depends on the data being modeled. A commonly used kernel is the radial basis function kernel (RBF):

$$k(x, x') = \sigma_f^2 \exp \left[ - \frac{(x - x')^2}{2l^2} \right]. \quad (1)$$

It is characterized by two hyperparameters,<fnref n="16" /> $\sigma_f$ and $l$, where $\sigma_f$ describes the vertical extent of the function and the length scale $l$ describes how quickly the correlation between two points reduces as their distance increases (larger $l$ implies a smoother function). The hyperparameters of the model were fine-tuned within the scikit-learn module, which minimizes the log-likelihood on the available data. The best-fit GPR model can be used to generate predictions $f$ and their associated uncertainty var($f$) as follows:

$$f = K_*^T (K + \sigma_n^2 I)^{-1} y, \quad (2)$$

$$\operatorname{var}(f) = \operatorname{diag}[K_{**} - K_{*}^{T}(K + \sigma_{n}^{2}I)^{-1}K_{*}], \quad (3)$$

where $I$ is the identity matrix, $y$ are the observed data points (in our case the timing residuals), the elements of covariance matrices are $K = k(x, x)$, $K_* = k(x, x')$, $K_{**} = k(x', x')$, and the apex “T” indicates the transpose matrix. For more details on the GPR method and its applications in pulsar timing we refer

the readers to C. E. Rasmussen & C. K. I. Williams (2006), P. R. Brook et al. (2016), and K. M. Rajwade et al. (2022). In our case, the timing residuals were best modeled with two RBF kernels in addition to a white-noise kernel that accounts for uncertainties $\sigma_n$ in the timing residuals. The optimized hyperparameters were $\sigma_{f_1} = 0.254 \pm 0.01$ cycles, $l_1 = 11.5 \pm 0.3$ days, $\sigma_{f_2} = 26.5_{-2}^{+1}$ cycles, $l_2 = 37.6_{-1}^{+0.4}$ days, and $\sigma_n = (1.7 \pm 0.4) \times 10^{-2}$ cycles.

We remark that, although the optimized GPR model does not provide any physical insights, it provides a means to evaluate derivatives of the timing residuals, which can be used to compute the observable parameters $\nu$ and $\dot{\nu}$. The derivatives and their uncertainties were evaluated analytically using the formalism outlined by P. R. Brook et al. (2016). The derived evolution in $\nu$ and $\dot{\nu}$ are indicated in the middle and bottom panels of Figure 1 and in Figure 2, suggesting that these values are generally consistent with those derived using the phase-coherent timing approach in the previous section and by C.-P. Hu et al. (2023).

Our analysis suggests that the onset of the X-ray outburst (in the period between 2020 May and June) triggered a dramatic change in the spin parameters, especially with the magnitude of $\dot{\nu}$ showing a relative increase (i.e., $\frac{\Delta\dot{\nu}}{\dot{\nu}}$) of more than 7%, before reaching values close to the pre-outburst level a few months later. Possible physical explanation for such a complex variation in the spin parameters will be discussed further in Section 7.

## 3.3. Pulse Profile Analysis

Obtaining local measurements of the spin frequency and spin-down rate via GPR and TEMPO (shown in Figure 1) enabled us to compute the pulse profiles by folding on the optimal ephemeris using routines in Stingray (D. Huppenkothen et al. 2019).<fnref n="17" /> These profiles are shown in Figures 3 and 4 for the four separate instruments. In the case of NICER and Swift the pulse profiles from several observations (indicated at the top of each figure) were combined to produce a stacked profile. The pulse profiles generally appear sinusoidal with minimal contribution from higher-order harmonics >4. However, we note that for the period just after the 2020 outburst (i.e., MJD 59,055–59,085), the stacked NICER profile features a shoulder to the main sinusoidal component, which was originally found by C.-P. Hu et al. (2023).

We determined the long-term evolution of the pulsed flux as illustrated in Figure 1 by performing phase-resolved spectroscopy, which is discussed further in Section 5. It is apparent that the pulsed flux increases notably on 2020 June 26 (see also L. Kuiper et al. 2020 and C.-P. Hu et al. 2023), which is coincident with a (spin-up) timing glitch and a millisecond X-ray burst that was observed a few days later via Swift-BAT (i.e., on MJD 59,063; L. Kuiper et al. 2020).

We also examined how the pulsed flux of the source varies as a function of energy for the XMM-Newton observations, as indicated in Figure 5. It should be noted that we compute the pulsed flux as a function of energy rather than the pulsed fraction, since the total flux (in this case) does not discount the contribution from the PWN, thus biasing its estimates. That said, the pulsed flux was estimated using phase-resolved spectroscopy (see Section 5) and hence subtracts the emission

<footnote n="16">¹⁶ They are defined as hyperparameters, since they are only used in training the regression model, rather than describing any physically relevant properties of the timing residuals.</footnote>

<footnote n="17">¹⁷ https://docs.stingray.science/index.html</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top panel: evolution of the pulsed flux (in the 1–12 keV energy range) covering the periods before (left blue-shaded region), during (red-shaded region), and after (right green-shaded region) the 2020 August outburst of PSR J1846−0258. The solid black line denotes the epoch when the Swift Burst Alert Telescope (BAT) and Fermi Gamma-ray Burst Monitor (GBM) both confirmed the detection of a short (<100 ms) burst from the source, while the dashed black line indicates a slightly earlier period during which a spin-up glitch/irregularity was found to occur with respect to the pre-outburst spin evolution. The pulsed flux was computed using phase-resolved spectroscopy after deriving the spin parameters $\nu$ and $\dot{\nu}$ via the tempo-derived phase-coherent timing approach (see Section 3.1) and Gaussian process regression (GPR; see Section 3.2). Middle panel: evolution in the spin frequency relative to the pre-outburst timing solution (i.e., before MJD 59,026) $\Delta\nu$ derived from computing the derivative of the best-fit Gaussian process model to the timing residuals. Bottom panel: evolution in the frequency derivative $\dot{\nu}$. For all three panels, the errors are quoted at the $1\sigma$ confidence level.</figcaption></figure>

from the PWN. The pulsed flux appears to show curvature at energies below 6 keV, while increasing to higher values above 10 keV. This clearly implies the presence of two spectral components to the pulsed emission (i.e., a blackbody and a power law; see Section 5 for more details on the modeling of this emission).

Finally, we illustrate how the morphology of the XMM-Newton pulse profiles vary with energy in Figure 4. As shown in the figure, we do not observe a significant change in the morphology of the pulse profiles or the occurrence of any phase shifts of the primary sinusoidal component with increasing energy. For some observations the pulse profile does appear to be more complex than a sinusoid, but no clear trends are observed as a function of energy or observation date.

## 4. Diffuse Emission

To analyze the extended emission surrounding the source (see Figure 6), we extracted surface brightness profiles of the three XMM-Newton observations in three different energy bands (0.3–2 keV, 2–5 keV, and 5–8 keV) out to a radius of 250" (in increments of 2"). These are shown in Figure 7 and were fitted with an analytical King point-spread function model plus a constant background term $b_0$ such that the surface

brightness $S_b$ was modeled as

$$S_b = \frac{a_0}{(1 + (r/r_c)^2)^\alpha} + b_0, \quad (4)$$

where $r$ is the distance from the center of the source, and $a_0$, $r_c$, $\alpha$, and $b_0$ are free parameters determined by the fit. While fitting the radial profiles with this model, we excluded the radial region between 35" and 160" in order to avoid the contribution from the Kes 75 SNR. However, for visual purposes we overlay the emission from the remnant in Figure 7 to illustrate the clear excess in $S_b$ in the abovementioned radial regions. The excess appears very strong in the 0.3–2 keV and 2–5 keV energy bands, but decreases in strength above 5 keV, illustrating the energy dependence of the emission from the remnant. As expected, there appears to be little variability in the excess emission between the three observations for all energy bands.

## 5. Phase-resolved X-Ray Spectroscopy

We made use of the NICER, NuSTAR, and XMM-Newton EPIC-pn observations to characterize the spectral properties of the pulsed emission over the energy energy range of 1–70 keV using the Xspec software package (K. A. Arnaud 1996). In order to isolate the pulsed spectral component, we used the spin
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison between the values of the spin frequency $\Delta\nu$ (relative to the pre-outburst timing model; left) and the spin-down rate $\dot{\nu}$ (right) obtained in our study through phase-coherent timing (black points) and Gaussian process regression (GPR). We also show the values reported by C.-P. Hu et al. (2023, solid blue line).</figcaption></figure>

(The data used to create this figure are available in the online article.)

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Stacked NICER (blue), NuSTAR (magenta), and Swift (red) pulse profiles for the periods before (MJD 58,352–58,996), during (MJD 59,055–59,085, MJD 59,087–59,122, and MJD 59,125–59,171) and after (MJD 59,281–59,435) the 2020 outburst of PSR J1846–0258. The pulse profiles were combined over the observational epochs indicated at the top of each profile in the energy ranges of 1–12 keV (for NICER), 3–79 keV (for NuSTAR), and 1–10 keV for Swift.</figcaption></figure>

parameters derived from the timing analysis discussed in the previous section to perform phase-resolved spectroscopy. For the NICER observations, we implemented the *photonphase*

tool to assign a pulse phase to each event for a given observation based on the available timing ephemeris, and extracted the spectra for those intervals corresponding to the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Energy-resolved pulse profiles from the three post-outburst XMM-Newton observations, with the red solid line displaying the best-fit sinusoidal model (with at most two harmonic components). The number of phase bins for each profile was chosen based on the total number of counts available. Observational dates are indicated in the top panel of each plot.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Pulsed flux as a function of energy for the three post-outburst XMM-Newton observations.</figcaption></figure>

brightest (on) and faintest (off) quarters of the pulse cycle using *niextract-events* and *niextspect*. This was followed by the generation of the ancillary response files (ARFs) and response matrix files (RMFs) via the *nicerarf* and *nicerrmf* tools. We note that care was taken to ensure that the nonbarycentered event files were used to generate the ARFs and RMFs, since it is known that the NICER pipeline could fail to produce the correct responses for barycentered event files. The pulse phase column generated on the barycentered event files were copied over to the nonbarycentered files using

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The XMM-Newton EPIC-pn image of the field containing the high-B pulsar PSR J1846−0258, its PWN, and the spatially resolved SNR Kes 75. The circular black solid line indicates the extraction region used for the phase-resolved spectral analysis of the source.</figcaption></figure>

`ftpaste`, and appropriate phase filters were then applied to extract the pulse-on and pulse-off spectra. The `nicerarf` tool requires the source position as input, since the detector quantum efficiency depends on the off-axis angle, and so we used the source position available from the Chandra observations. During the process, we also used the filter file for each observation containing key information regarding various screening criteria including passages of the spacecraft through the SAA. The pulsed spectra were finally extracted by subtracting the pulse-off spectrum from the pulse-on spectrum, and were grouped using an optimal binning algorithm from J. S. Kaastra & J. A. M. Bleeker (2016) via the `ftgrouppha` tool, with at least 200 counts per energy bin. We further placed
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Surface brightness (SB) profiles of the three XMM-Newton observations in three different energy bands up to a radial distance of 250" from the position of PSR J1846−0258. In each case, the radial profile was fit with a King point-spread function model (after excluding the radial range of 35"−160" to avoid the contribution from the SNR; see text; blue lines) with the addition of a constant background (black solid lines) in order to analyze the extension of the source. The residuals to the best fit are shown in the bottom panel.</figcaption></figure>

the pulsed NICER spectra into eight distinct groups, chosen to correspond to eight different flux levels spanning the entire duration of the outburst. Spectra in each group were combined using the tool `addspec`.

In order to process the NuSTAR data for phase-resolved spectroscopy, we generated Good Time Interval (GTI) files corresponding to the brightest and faintest quarters of the pulse cycle. The pulsed spectra were subsequently extracted with the `nuproducts` tool, adopting a source region size as discussed in Section 2. The spectra were rebinned using `ftgrouppha`, with at least 50 counts per energy bin.

For XMM-Newton, the spectral extraction procedure was broadly similar. We selected source and background events according to the procedure discussed in Section 2. The `photonphase` tool was used to assign a pulse phase to each event by folding the data on the best-fit ephemeris. We then selected the quarters of the pulse cycle corresponding to the brightest and faintest emission to generate a “pulse-on” and a “pulse-off” spectrum, and subtracted the latter from the former. We generated ARFs and RMFs using the `arfgen` and `rmfgen` tasks, respectively. All spectra were rebinned using the `ftgrouppha` tool adopting at least 200 counts per energy bin, and the spectra from the three instruments were read into `Xspec`. Finally, we found that the Swift data were inappropriate for phase-resolved analysis due to the paucity of counts, although they were used to estimate the pulsed flux.

In Figure 8, we present the pulsed spectra over the 1–70 keV energy range. We proceeded to fit the spectra collectively with an absorbed single-temperature blackbody plus power-law

model. The blackbody temperature and normalization were allowed to vary for each spectrum, while the column density was tied between the different observations. The power-law slope and normalization were allowed to vary only for the NuSTAR spectra, but for the XMM-Newton and NICER spectra, due to their limited energy range, the power-law parameters were fixed to (weighted) average values found for the four NuSTAR spectra. It should be noted, however, that the power-law component remains largely unchanged, as illustrated in Figure 8. In all cases, we used the abundances by J. Wilms et al. (2000). The best-fit parameters of the spectral analysis are reported in Table 1, with an overall goodness-of-fit that is statistically acceptable (with $\chi^2$/dof = 1372.0/1081).

## 6. Radio Observations

PSR J1846−0258 was observed in radio with the Parkes-Murriyang telescope on 2020 August 5. Observations, obtained through a target of opportunity observation (ID PX065), were performed with the Ultra-Wide-bandwidth, Low-frequency receiver (UWL; G. Hobbs et al. 2020), having a total bandwidth of 3328 MHz centered at 2368 MHz. The band was split into 3328 frequency channels. Full Stokes data were recorded and 8 bit sampled every 512 µs for 3.9 hr, starting at UT 08:37. The data were searched for periodic pulsation and for single dispersed pulses.

The search for a periodic signal has been performed both in the entire observational bandwidth and in three separate subbands (702–1000 MHz, 1000–2000 MHz, and 2000–4032 MHz).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top panel: absorbed ($N_H = (5.5 \pm 0.3) \times 10^{22} \text{ cm}^{-2}$) unfolded (deconvolved) pulsed energy spectra of PSR J1846−0258 in the aftermath of the 2020 outburst using NICER, NuSTAR, and XMM-Newton observations (see text). For display purposes, the NICER spectra were separated into seven groups, where in each group the spectra were stacked using the *addspec* tool (for the observation date ranges indicated in the legend). Similarly, the three XMM-Newton spectra were also combined. However, for spectral fitting we used the original spectra for each observation (see Table 1). Each spectrum was fit with an absorbed blackbody plus power-law model with the column density tied between the different observations. The results clearly demonstrate the gradual fading of the thermal component (below 10 keV) during the outburst, while the power-law component (above 8 keV) remains unchanged. The detailed best-fit parameters are reported in Table 1. The residuals of the best-fit model are also shown in the bottom panel of the figure.</figcaption></figure>

Data were folded with the partially coherent solution covering the radio epoch and with dispersion measures (DMs) from 0 to $2000 \text{ pc cm}^{-3}$. Folded archives with 1 MHz frequency resolution and 512 time bins across the folded profile were created with a DM step of 200 pc cm$^{-3}$. After radio frequency interference (RFI) removal, each archive was searched over a period range spanning $\pm 3.5 \mu\text{s}$ and $\pm 100 \text{ pc cm}^{-3}$ around the nominal values.

The single-pulse search was done by implementing an iterative subbanded search (similarly to P. Kumar et al. 2021), in addition to the search in the full observational band. The reason for subbanding is that repeating FRBs have been observed as band-limited sources, meaning that the emission is typically limited to as little as ~100 MHz, especially at the lowest frequencies (P. Kumar et al. 2021), and this would result in a large loss of sensitivity (or of the burst altogether) in the full-band search. In the subbanded search, we searched for bursts by inspecting iteratively smaller bandwidths of widths [2 × 1664, 4 × 832, 8 × 416, 16 × 208, 32 × 104] MHz. For each subband, we also considered overlapping adjacent subbands (by shifting the bands by half widths) in order not to miss events between two adjoining subbands. This results to a total of [1 × 3328, 3 × 1664, 7 × 832, 15 × 416, 31 × 208, 63 × 104] MHz bands processed. Data were dedispersed from 0 to $2000 \text{ pc cm}^{-3}$ and we searched for excesses by convolving

the frequency-collapsed time series with a top-hat function with trial widths logarithmically spaced and ranged up to a maximum width of 512 bins.

No clear sign of either periodic radio emission or sporadic single pulses was found above a signal-to-noise ratio of 8 and 10, respectively (see also the works by H. Blumer et al. 2021 at 2 GHz, M. B. Mickaliger et al. 2020 at 1.5 GHz, and W. A. Majid et al. 2020 at 8.3 GHz and 13.9 GHz). From the single-pulse search we redetected the RRAT J1846-0257 (M. A. McLaughlin et al. 2009), whose position is ~1.3 from the nominal center of the telescope beam pointing, with an FWHM beamwidth of ~4.9 at 4032 MHz. By exploiting the standard radiometer equation for folded pulsar profiles and single pulses (J. M. Cordes & M. A. McLaughlin 2003; D. Lorimer & M. Kramer 2004), we can set our upper limits both for radio periodic emission and single bursts. We assume for the UWL a telescope gain of $G = 1.8 \text{ Jy K}^{-1}$ and an average system temperature without the sky contribution of ~22 K (see Figure 4 from G. Hobbs et al. 2020). The frequency-dependent sky temperatures at the source pointing were computed assuming a pivotal temperature of 342.3 K at 408 MHz (C. G. T. Haslam et al. 1982) and scaled with a spectral index of -2.7.
<table><caption>**Table 1**
The Best-fit Parameters Resulting from Fitting the Pulsed Spectra of All Three Instruments</caption><thead><tr><th>Instrument</th><th>Epoch (MJD)</th><th>kT<sub>bb</sub> (keV)</th><th>R<sub>bb</sub><sup>a</sup> (km)</th><th>Total Flux<sup>b</sup> (10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>BB Flux<sup>b</sup> (10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>χ<sup>2</sup>/dof<sup>c</sup></th></tr></thead><tbody><tr><td>NICER</td><td>59,026</td><td>1.0 ± 0.3</td><td>0.5 ± 0.3</td><td>1.3 ± 0.1</td><td>0.6 ± 0.1</td><td>43.4/44</td></tr><tr><td>NICER</td><td>59,055</td><td>1.1 ± 0.2</td><td>1.0<sup>+0.5</sup><sub>-0.4</sub></td><td>3.6 ± 0.3</td><td>3.1 ± 0.5</td><td>18.2/18</td></tr><tr><td>NICER</td><td>59,056</td><td>1.1 ± 0.2</td><td>0.9<sup>+0.5</sup><sub>-0.4</sub></td><td>3.4 ± 0.2</td><td>2.8 ± 0.4</td><td>30.3/19</td></tr><tr><td>NICER</td><td>59,064</td><td>0.95 ± 0.1</td><td>1.0 ± 0.4</td><td>2.6 ± 0.1</td><td>2.0 ± 0.2</td><td>64.1/38</td></tr><tr><td>NICER</td><td>59,067</td><td>1.0 ± 0.1</td><td>0.8 ± 0.3</td><td>2.5 ± 0.1</td><td>1.9 ± 0.3</td><td>49.9/30</td></tr><tr><td>NICER</td><td>59,068</td><td>1.0 ± 0.2</td><td>0.8<sup>+0.6</sup><sub>-0.4</sub></td><td>2.3 ± 0.2</td><td>1.7 ± 0.4</td><td>18.5/20</td></tr><tr><td>NICER</td><td>59,070</td><td>1.0 ± 0.1</td><td>0.8 ± 0.3</td><td>2.4 ± 0.1</td><td>1.7 ± 0.3</td><td>37.9/41</td></tr><tr><td>NICER</td><td>59,072</td><td>1.0 ± 0.2</td><td>0.8<sup>+0.5</sup><sub>-0.4</sub></td><td>2.1 ± 0.2</td><td>1.4 ± 0.4</td><td>16.5/24</td></tr><tr><td>NICER</td><td>59,073</td><td>0.9 ± 0.2</td><td>0.9<sup>+0.4</sup><sub>-0.3</sub></td><td>2.1 ± 0.1</td><td>1.5 ± 0.2</td><td>42.2/33</td></tr><tr><td>NICER</td><td>59,075</td><td>1.0 ± 0.2</td><td>0.7<sup>+0.5</sup><sub>-0.3</sub></td><td>2.0 ± 0.1</td><td>1.4 ± 0.3</td><td>24.5/22</td></tr><tr><td>NICER</td><td>59,079</td><td>0.9 ± 0.2</td><td>1.0<sup>+0.5</sup><sub>-0.4</sub></td><td>2.4 ± 0.2</td><td>1.7 ± 0.3</td><td>36.6/28</td></tr><tr><td>NICER</td><td>59,081</td><td>1.1 ± 0.2</td><td>0.7 ± 0.3</td><td>3.0 ± 0.1</td><td>2.3 ± 0.4</td><td>44.0/33</td></tr><tr><td>NICER</td><td>59,086</td><td>1.0 ± 0.2</td><td>0.8<sup>+0.4</sup><sub>-0.3</sub></td><td>2.3 ± 0.5</td><td>1.7 ± 0.4</td><td>38.5/35</td></tr><tr><td>NICER</td><td>59,088</td><td>0.9 ± 0.2</td><td>0.8<sup>+0.4</sup><sub>-0.3</sub></td><td>2.0 ± 0.1</td><td>1.3 ± 0.3</td><td>86.2/45</td></tr><tr><td>NICER</td><td>59,092</td><td>0.8 ± 0.2</td><td>0.9<sup>+0.7</sup><sub>-0.5</sub></td><td>1.6 ± 0.1</td><td>0.9 ± 0.3</td><td>68.1/41</td></tr><tr><td>NICER</td><td>59,098</td><td>1.0 ± 0.2</td><td>0.7<sup>+0.5</sup><sub>-0.3</sub></td><td>1.8 ± 0.1</td><td>1.2 ± 0.2</td><td>65.0/48</td></tr><tr><td>NICER</td><td>59,102</td><td>0.7 ± 0.2</td><td>1.2<sup>+1</sup><sub>-0.8</sub></td><td>1.4 ± 0.2</td><td>0.8 ± 0.3</td><td>46.9/40</td></tr><tr><td>NICER</td><td>59,122</td><td>1.2 ± 0.4</td><td>0.35<sup>+0.5</sup><sub>-0.2</sub></td><td>1.2 ± 0.1</td><td>0.6 ± 0.3</td><td>49.8/41</td></tr><tr><td>NICER</td><td>59,137</td><td>0.9 ± 0.2</td><td>0.5<sup>+0.6</sup><sub>-0.3</sub></td><td>1.2 ± 0.1</td><td>0.5 ± 0.3</td><td>52.5/28</td></tr><tr><td>NICER</td><td>59,145</td><td>0.6 ± 0.4</td><td>0.8<sup>+6</sup><sub>-0.6</sub></td><td>0.9 ± 0.2</td><td>&lt;0.7</td><td>66.3/60</td></tr><tr><td>NICER</td><td>59,154</td><td>0.2 ± 0.1</td><td>&lt;170</td><td>1.5 ± 0.2</td><td>&lt;1.5</td><td>92.3/49</td></tr><tr><td>NICER</td><td>59,165</td><td>0.2 ± 0.1</td><td>&lt;3670</td><td>1.5 ± 0.2</td><td>&lt;1.6</td><td>53.2/40</td></tr><tr><td>NICER</td><td>59,174</td><td>0.4<sup>+0.1</sup><sub>-0.3</sub></td><td>&lt;18000</td><td>0.6 ± 0.2</td><td>&lt;0.2</td><td>28.8/16</td></tr><tr><td>XMM-Newton</td><td>59,109</td><td>&lt;0.9</td><td>&lt;0.5</td><td>0.9 ± 0.2</td><td>&lt;0.2</td><td>75.6/74</td></tr><tr><td>XMM-Newton</td><td>59,125</td><td>0.7<sup>+0.3</sup><sub>-0.2</sub></td><td>0.5<sup>+1</sup><sub>-0.3</sub></td><td>1.0 ± 0.1</td><td>0.3 ± 0.1</td><td>77.9/73</td></tr><tr><td>XMM-Newton</td><td>59,145</td><td>0.9<sup>+0.6</sup><sub>-0.4</sub></td><td>0.3<sup>+0.6</sup><sub>-0.1</sub></td><td>0.9 ± 0.1</td><td>0.2 ± 0.1</td><td>84.0/84</td></tr></tbody></table><table><thead><tr><th>Instrument</th><th>Epoch</th><th>kT<sub>bb</sub> (keV)</th><th>R<sub>bb</sub> (km)</th><th>Total flux (10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Photon index</th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td>NuSTAR</td><td>590,66</td><td>0.9 ± 0.1</td><td>0.8<sup>+1</sup><sub>-0.4</sub></td><td>3.1 ± 0.4</td><td>1.0 ± 0.7</td><td>22.1/18</td></tr><tr><td>NuSTAR</td><td>59,082</td><td>0.9 ± 0.1</td><td>1.0<sup>+0.5</sup><sub>-0.3</sub></td><td>3.5 ± 0.2</td><td>&lt;1.5</td><td>26.7/26</td></tr><tr><td>NuSTAR</td><td>59,109</td><td>...</td><td>...</td><td>3.5 ± 0.1</td><td>1.4 ± 0.2</td><td>26.0/22</td></tr><tr><td>NuSTAR</td><td>59,132</td><td>...</td><td>...</td><td>2.5 ± 0.3</td><td>1.3 ± 0.2</td><td>32.2/35</td></tr></tbody><tfoot><tr><td colspan="6"></td><td>1372.0/1081</td></tr></tfoot></table>

Notes. The column density was tied between the different observations, yielding $N_H = (5.5 \pm 0.3) \times 10^{22} \text{ cm}^{-2}$. All errors are quoted at the $1\sigma$ confidence level.

a Assuming a source distance of 5.8 kpc (J. P. W. Verbiest et al. 2012).

b Energy range: 1–10 keV.

c dof: degrees of freedom.

d Energy range: 3–79 keV.

For the putative periodic radio emission, assuming a fiducial period of $P = 0.326$ s and a duty cycle of 10%, we report the peak flux density upper limits in the four bands processed in Table 2.

Table 2 also reports our fluence density upper limits for the emission of single pulses in the full band and the subbands (the upper limits are quoted at the $1\sigma$ confidence level). Here, we assumed a pulse/burst with 1 ms duration and considered, for each chosen bandwidth, the lowest frequency subband (as this is the one with the highest sky temperature). The limits in the subbanded search are, as expected, higher than the limit in the full band. This limit, however, works under the assumption that the burst is visible within the full bandwidth and that it has a flat spectrum. In the hypothesis that the burst is band limited, the subbanded search can become more sensitive as the average signal-to-noise ratio will be less affected by the noise contribution (for details, see M. Trudu et al. 2024, in preparation).

## 7. Discussion

The source PSR J1846−0258 is one of the best examples of a highly energetic pulsar showing both rotation and magnetic-powered emission during its bursting period. This pulsar, together with PSR J1119−6127 and the low-field magnetars, suggests how different pulsar classes can be unified based on their age, magnetic field strength, and geometric configuration at birth. In this work, we studied in detail the 2020 outburst of PSR J1846−0258. Among other findings, we found indications for an intriguing oscillation of the spin-frequency derivative in the immediate aftermath of the 2020 outburst (see Figures 1 and 2). This is rather unexpected as quasi-periodic oscillations in the timing residuals of isolated pulsars are the result of timing noise that occurs on timescales of several years (see, e.g., G. Hobbs 2009). However, in this case, the timescales are significantly shorter.

Considering a time variable coupling of the superfluid NS core and crust (i.e., the model by E. Güercinoğlu et al. 2023),
<table><caption>**Table 2**
Radio Emission Upper Limits (Quoted at the 1σ Confidence Level)</caption><thead><tr><th>Periodic Emission</th><th>Frequency Band (MHz)</th><th>S<sub>ν</sub> (mJy)</th></tr></thead><tbody><tr><td rowspan="4"></td><td>702–4032</td><td>0.01</td></tr><tr><td>702–1000</td><td>0.03</td></tr><tr><td>1000–2000</td><td>0.02</td></tr><tr><td>2000–4032</td><td>0.01</td></tr><tr><th>Single pulses</th><th>Frequency band (MHz)</th><th>F<sub>ν</sub> (Jy ms)</th></tr><tr><td></td><td>702–4032</td><td>0.2</td></tr><tr><td></td><td>702–2366</td><td>0.3</td></tr><tr><td></td><td>702–1534</td><td>0.9</td></tr><tr><td></td><td>702–1118</td><td>2</td></tr><tr><td></td><td>702–954</td><td>3</td></tr><tr><td></td><td>702–806</td><td>4</td></tr></tbody></table>

Note. For the periodic emission we report the peak flux density $S_ν$ upper limit, whereas for single pulses we report the fluence density $F_ν$ (see Section 6 for more details).

we attempt to constrain the timescale of the quasi-periodicity by modeling the evolution of $\dot{\nu} = \dot{\Omega}/2\pi$ with a combination of a damped simple harmonic oscillator (SHO) and regular spin-down, in line with the following expression:

$$ \begin{aligned} \dot{\Omega}(t) = & 2\pi\dot{\nu}_0 + 2\pi\ddot{\nu}_0 t \\ & + A\Omega_0 \exp\left(-\frac{t}{2\tau}\right) \cos[\Omega_0 t + \phi] \\ & - \frac{A}{2\tau} \exp\left(-\frac{t}{2\tau}\right) \sin[\Omega_0 t + \phi] \quad \text{with} \\ \Omega_0 = & \omega_0 \left[1 - \left(\frac{1}{2\tau\omega_0}\right)^2\right]^{1/2}, \end{aligned} \quad (5) $$

where the symbols are defined in Table 3 and have the same meanings as in E. Güercinoğlu et al. (2023).

We also consider the precession model of G. Ashton et al. (2017; triggered plausibly due to magnetic readjustments in the NS crust at the outburst onset, leading to deformation of the crust) in order to fit the periodic trend in the spin-down rate. We allow the precession angle to vary secularly with time (see Equations (2)–(6) and (10) of G. Ashton et al. 2017 for full details on the model).

The above two models were fit to the data using least-squares minimization via the `kmpfit` package (J. P. Terlouw & M. G. R. Vogelaar 2014), but the resulting fits (after convergence) were found to be statistically unacceptable. This is apparent from the residuals to the best-fit models that clearly show considerable scatter (see the bottom panel of Figure 9). Moreover, performing a runs test rejects the two models at a significance of above $5\sigma$ (with p-values below $10^{-7}$). In order to improve the goodness-of-fit and extract parameter uncertainties, we fit the models to a smaller section of the data (i.e., MJD 59,055–59,135). We show the corresponding best-fit parameters in Tables 3 and 4. We note that, while the above models do not formally describe the full data set, they provide a crude estimate of the oscillation period $P_{\text{osc}} = \frac{2\pi}{\Omega_0} \sim 50–60$ days and the damping timescale $\tau \sim 20$ days (in the context of the damped harmonic oscillator) and the ellipticity $\epsilon$ of the NS (for the precession model), which we discuss below.

<table><caption>**Table 3**
The Best-fit Parameters for Modeling the Evolution of $\dot{\nu}$ with a Damped Harmonic Oscillator</caption><thead><tr><th>Parameter</th><th>Best-fit Value</th></tr></thead><tbody><tr><td>Validity period (MJD)</td><td>59,055–59,135</td></tr><tr><td>φ (reference phase)</td><td>-0.7 ± 0.1</td></tr><tr><td>ω<sub>0</sub></td><td>0.120 ± 0.002</td></tr><tr><td>Damping timescale, τ (days)</td><td>20 ± 2</td></tr><tr><td>Amplitude, A (rad s<sup>-1</sup>)</td><td>(2.7 ± 0.2) × 10<sup>-6</sup></td></tr><tr><td>ν̇<sub>0</sub> (Hz s<sup>-1</sup>)</td><td>- (7.06 ± 0.01) × 10<sup>-11</sup></td></tr><tr><td>ν̇<sub>0</sub> (Hz s<sup>-2</sup>)</td><td>(1.9 ± 0.2) × 10<sup>-19</sup></td></tr><tr><td>Oscillation frequency, Ω<sub>0</sub> (rad s<sup>-1</sup>)</td><td>(1.35 ± 0.03) × 10<sup>-6</sup></td></tr><tr><td>χ<sup>2</sup>/dof</td><td>51.5/46</td></tr></tbody></table>

We suggest leaving a detailed modeling of the spin-down rate to a future work, while noting the main points below.

In the context of the damped SHO model, the oscillation period of ~50–60 days could be consistent with Ekman oscillations (due to nonspherical geometry of the superfluid vortex lines) with $T_{\text{Ekman}} = \frac{R}{(\kappa\Omega^2)^{1/2}}$, where R is the NS radius, $\Omega$ is the angular rotational velocity, and $\kappa = 2 \times 10^{-3} \text{ cm}^2 \text{ s}^{-1}$ is the quantized vorticity attached to each vortex line (see Section 3 of E. Güercinoğlu et al. 2023). Moreover, the damping timescale $\tau$ of the spin-down rate is related to the linear creep timescale of the superfluid vortices, which depends on the microphysical properties of the NS crust (see Equation (5) of E. Güercinoğlu et al. 2023). The observed $\tau \sim 20$ days seems plausible if the surface temperature of the NS is ~0.015 keV. We note that this is much smaller than the blackbody temperature observed during late stages of the outburst in this work (although a more precise measurement of the surface temperature would require optical and UV observations and precise constraints on the hydrogen column density). In computing the expected damping timescale, we assume typical parameter values for other microphysical quantities of the NS (e.g., the age-dependent critical angular velocity between normal matter and superfluid crust $\omega_{\text{cr}} \sim 5 \times 10^{-3} \text{ rad s}^{-1}$, the pinning energy between superfluid vortices and lattice nuclei $E_{\text{P}} \sim 10$ keV, and microscopic vortex velocity around nuclei $v_0 \sim 10^7 \text{ cm s}^{-1}$; see M. A. Alpar et al. 1989; B. Haskell & A. Melatos 2015).

In the context of the precession model, the magnetic ellipticity of the NS is given by $\epsilon_B = \frac{B^2 R^4}{GM^2}$ and is about $5.7 \times 10^{-8}$ for this source assuming that the surface dipolar magnetic field $B$ (~$5 \times 10^{13}$ G) can be estimated in the standard manner via the observed spin period and spin-down rate at outburst onset for a canonical mass of $1.4 M_\odot$ and radius $R$ of 10–12 km. The estimated ellipticity seems consistent with the ellipticity inferred by modeling the spin-down rate variability (see Table 4). However, we note that the $\dot{\nu}$ modulation does not appear to be correlated with (periodic) changes in the source flux, which could disfavor an interpretation in terms of precession (although this statement may depend on the size of the precession angle).

Furthermore, the magnetospheric twist model developed by A. M. Beloborodov (2009) is often invoked to explain the spin-down and radiative properties of a large number of magnetars. In this model, the crust of the NS is prone to azimuthal displacements resulting from star quakes, which is expected to implant twists in the current-carrying bundle of poloidal field
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The evolution in the frequency derivative $\dot{\nu}$ modeled with the combination of spin-down terms and damped harmonic oscillator (blue curve; E. Güercinoğlu et al. 2023) and free precession (red curve; G. Ashton et al. 2017), as indicated in the top panel. The residuals (data – model) are shown in the bottom panel, which clearly demonstrate significant scatter (indicating that these models are likely unsuitable to describe the full complex $\dot{\nu}$ evolution). To obtain a statistically acceptable fit, we fit the models to a smaller section of the data (i.e., MJD 59,055–59,135). We find for the damped SHO, $\chi^2/\text{dof} = 51.5/46$, while for the free precession model, $\chi^2/\text{dof} = 50.5/45$.</figcaption></figure>

<table><caption>**Table 4**
The Best-fit Parameters for Modeling the Evolution of $\dot{\nu}$ with Free Precession</caption><thead><tr><th>Parameter</th><th>Best-fit Value</th></tr></thead><tbody><tr><td>Validity period (MJD)</td><td>59,055–59,135</td></tr><tr><td>T<sub>ref</sub> (days after outburst)</td><td>59,076.9 ± 0.3</td></tr><tr><td>Wobble angle θ<sub>0</sub> (rad)</td><td>(7.3 ± 0.1) × 10<sup>-5</sup></td></tr><tr><td>θ̇<sub>0</sub> (rad s<sup>-1</sup>)</td><td>-(2.0 ± 0.1) × 10<sup>-6</sup></td></tr><tr><td>χ (radians)</td><td>9.4285 ± 0.00001</td></tr><tr><td>ε</td><td>(6.0 ± 0.2) × 10<sup>-8</sup></td></tr><tr><td>ν̇<sub>0</sub> (Hz s<sup>-1</sup>)</td><td>-(7.048 ± 0.005) × 10<sup>-11</sup></td></tr><tr><td>ν<sub>0</sub> (Hz s<sup>-2</sup>)</td><td>(5.2 ± 0.3) × 10<sup>-19</sup></td></tr><tr><td>ψ<sub>0</sub></td><td>0</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>50.5/45</td></tr></tbody></table>

lines in the regions above the crust. The twist inflates the poloidal field lines, leading initially to a stronger dipole torque on the NS and thus increasing the spin-down rate following the outburst. However, these twists are gradually expected to decay due to the generation of screening currents (e.g., electron-positron pairs) in the magnetosphere, which oppose the currents maintaining the twist. This then implies that the spin-down torque acting on the NS decreases, and hence the spin-down rate also decreases. Qualitatively, the observed

evolution of $\dot{\nu}$ shown in Figure 1 could be consistent with the above model, although the presence of an oscillation in $\dot{\nu}$ during the outburst may be difficult to explain, and there is also the question of whether the timescale over which the $\dot{\nu}$ varies can be reasonably accounted for.

We also carried out radio observations of the source (along with the X-ray observations; see Table 5) with the Parkes-Murriyang telescope, setting upper limits for periodic radio emission and single pulses at frequencies centered at 2368 MHz, with a bandwidth of 3328 MHz. The radio-quiet nature of the source has previously been found at other frequencies by M. B. Mickaliger et al. (2020), H. Blumer et al. (2021), and W. A. Majid et al. (2020).

Moreover, phase-resolved X-ray spectroscopy carried out using all four instruments clearly demonstrates that the pulsed flux of the source (computed over the 1–12 keV energy range) begins to rise from the quiescent level at around MJD 59,026 (see Figure 1), peaks at MJD 59,056, and thereafter decreases exponentially before reaching pre-outburst levels 100 days later. Moreover, the source appears to undergo a second mini-outburst (C.-P. Hu et al. 2023) episode (lasting less than a day) at around MJD 59,082, when the decaying trend in the pulsed flux is interrupted by a sudden increase, as suggested by the NICER (Obs ID: 3033290116) and NuSTAR (Obs ID: 80602315004) observations.
<table><caption>**Table 5**
Observation Log</caption><thead><tr><th>Telescope</th><th>Obs ID</th><th>Start Day (YYYY-MM-DD)</th><th>Exposure<sup>a</sup> (ks)</th><th>Pulsed Flux<sup>b</sup> (10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>NICER</td><td>1033290117</td><td>2018-08-22</td><td>1.4</td><td>0.9 ± 0.6</td></tr><tr><td>NICER</td><td>1033290133</td><td>2018-11-10</td><td>5.9</td><td>0.4 ± 0.2</td></tr><tr><td>NICER</td><td>1033290135</td><td>2018-11-14</td><td>1.2</td><td>0.9 ± 0.4</td></tr><tr><td>NICER</td><td>1033290139</td><td>2018-11-18</td><td>1.3</td><td>0.6 ± 0.4</td></tr><tr><td>NICER</td><td>2516010101</td><td>2019-03-10</td><td>6.9</td><td>0.3 ± 0.2</td></tr><tr><td>NICER</td><td>2516010301</td><td>2019-03-13</td><td>4.9</td><td>0.3 ± 0.2</td></tr><tr><td>NICER</td><td>2516010401</td><td>2019-03-17</td><td>3.1</td><td>0.6 ± 0.3</td></tr><tr><td>NICER</td><td>2516010601</td><td>2019-05-06</td><td>4.4</td><td>0.5 ± 0.2</td></tr><tr><td>NICER</td><td>2516010801</td><td>2019-07-06</td><td>1.1</td><td>&lt;2</td></tr><tr><td>NICER</td><td>2516011102</td><td>2019-09-24</td><td>2.9</td><td>0.8 ± 0.5</td></tr><tr><td>NICER</td><td>2516011301</td><td>2019-09-28</td><td>2.6</td><td>0.9 ± 0.4</td></tr><tr><td>NICER</td><td>2516011701</td><td>2019-11-13</td><td>2.3</td><td>0.7 ± 0.4</td></tr><tr><td>NICER</td><td>2516011407</td><td>2019-02-14</td><td>2.8</td><td>0.9 ± 0.4</td></tr><tr><td>NICER</td><td>2516011406</td><td>2020-02-13</td><td>4.3</td><td>0.7 ± 0.3</td></tr><tr><td>NICER</td><td>2516011408</td><td>2020-02-15</td><td>6.3</td><td>0.4 ± 0.2</td></tr><tr><td>NICER</td><td>2516011502</td><td>2020-02-27</td><td>2.2</td><td>&lt;1</td></tr><tr><td>NICER</td><td>2516011503</td><td>2020-02-27</td><td>5.7</td><td>0.5 ± 0.3</td></tr><tr><td>NICER</td><td>3598010101</td><td>2020-03-26</td><td>6.1</td><td>0.5 ± 0.2</td></tr><tr><td>NICER</td><td>3598010302</td><td>2020-03-29</td><td>6.8</td><td>0.6 ± 0.2</td></tr><tr><td>NICER</td><td>3598010401</td><td>2020-04-02</td><td>2.6</td><td>0.6 ± 0.3</td></tr><tr><td>NICER</td><td>3598010402</td><td>2020-04-03</td><td>1.8</td><td>0.6 ± 0.4</td></tr><tr><td>NICER</td><td>3598010501</td><td>2020-04-29</td><td>9.4</td><td>0.6 ± 0.2</td></tr><tr><td>NICER</td><td>3598010601</td><td>2020-05-27</td><td>9.6</td><td>0.4 ± 0.2</td></tr><tr><td>NICER</td><td>3598010701</td><td>2020-06-26</td><td>7.3</td><td>1.6 ± 0.2</td></tr><tr><td>NICER</td><td>3598010801</td><td>2020-07-25</td><td>3.8</td><td>3.7 ± 0.3</td></tr><tr><td>NICER</td><td>3598010802</td><td>2020-07-25</td><td>4.3</td><td>3.6 ± 0.2</td></tr><tr><td>NICER</td><td>3033290101</td><td>2020-08-02</td><td>1.7</td><td>2.8 ± 0.4</td></tr><tr><td>NICER</td><td>3033290102</td><td>2020-08-02</td><td>9.0</td><td>2.5 ± 0.2</td></tr><tr><td>Swift</td><td>00032031217</td><td>2020-08-02</td><td>4.0</td><td>2.9 ± 0.8</td></tr><tr><td>NuSTAR</td><td>80602315002</td><td>2020-08-05</td><td>41.6</td><td>1.6 ± 0.2</td></tr><tr><td>NICER</td><td>3033290103</td><td>2020-08-05</td><td>5.2</td><td>2.5 ± 0.2</td></tr><tr><td>Swift</td><td>00032031218</td><td>2020-08-05</td><td>2.5</td><td>2.0 ± 0.9</td></tr><tr><td>NICER</td><td>3033290104</td><td>2020-08-06</td><td>5.4</td><td>2.2 ± 0.2</td></tr><tr><td>NICER</td><td>3033290105</td><td>2020-08-06</td><td>5.2</td><td>2.2 ± 0.2</td></tr><tr><td>NICER</td><td>3033290106</td><td>2020-08-08</td><td>4.7</td><td>2.3 ± 0.2</td></tr><tr><td>NICER</td><td>3033290107</td><td>2020-08-08</td><td>6.9</td><td>2.3 ± 0.2</td></tr><tr><td>NICER</td><td>3033290108</td><td>2020-08-10</td><td>6.0</td><td>2.1 ± 0.2</td></tr><tr><td>NICER</td><td>3033290109</td><td>2020-08-11</td><td>5.2</td><td>2.1 ± 0.2</td></tr><tr><td>NICER</td><td>3033290110</td><td>2020-08-12</td><td>4.5</td><td>2.2 ± 0.2</td></tr><tr><td>Swift</td><td>00032031219</td><td>2020-08-12</td><td>3.0</td><td>1.9 ± 0.8</td></tr><tr><td>NICER</td><td>3033290111</td><td>2020-08-13</td><td>3.3</td><td>2.1 ± 0.2</td></tr><tr><td>NICER</td><td>3033290112</td><td>2020-08-14</td><td>2.2</td><td>1.8 ± 0.3</td></tr><tr><td>NICER</td><td>3033290113</td><td>2020-08-17</td><td>3.5</td><td>2.4 ± 0.2</td></tr><tr><td>NICER</td><td>3033290114</td><td>2020-08-18</td><td>2.0</td><td>2.3 ± 0.3</td></tr><tr><td>NICER</td><td>3033290115</td><td>2020-08-19</td><td>3.9</td><td>2.0 ± 0.2</td></tr><tr><td>NuSTAR</td><td>80602315004</td><td>2020-08-20</td><td>56.2</td><td>2.5 ± 0.2</td></tr><tr><td>Swift</td><td>00032031221</td><td>2020-08-20</td><td>2.2</td><td>2.1 ± 1</td></tr><tr><td>NICER</td><td>3033290116</td><td>2020-08-21</td><td>3.1</td><td>4.2 ± 0.3</td></tr><tr><td>NICER</td><td>3598010901</td><td>2020-08-24</td><td>8.7</td><td>2.3 ± 0.2</td></tr><tr><td>NICER</td><td>3033290117</td><td>2020-08-26</td><td>4.7</td><td>2.0 ± 0.2</td></tr><tr><td>NICER</td><td>3033290118</td><td>2020-08-28</td><td>3.3</td><td>1.9 ± 0.3</td></tr><tr><td>NICER</td><td>3033290119</td><td>2020-08-30</td><td>3.8</td><td>1.5 ± 0.3</td></tr><tr><td>NICER</td><td>3033290120</td><td>2020-09-01</td><td>3.1</td><td>1.7 ± 0.3</td></tr><tr><td>NICER</td><td>3033290121</td><td>2020-09-03</td><td>1.3</td><td>1.4 ± 0.4</td></tr><tr><td>NICER</td><td>3033290122</td><td>2020-09-05</td><td>4.7</td><td>1.5 ± 0.2</td></tr><tr><td>NICER</td><td>3033290123</td><td>2020-09-07</td><td>4.5</td><td>2.2 ± 0.2</td></tr><tr><td>NICER</td><td>3033290124</td><td>2020-09-09</td><td>2.2</td><td>1.6 ± 0.3</td></tr><tr><td>NICER</td><td>3033290125</td><td>2020-09-11</td><td>2.2</td><td>1.5 ± 0.4</td></tr><tr><td>Swift</td><td>00032031225</td><td>2020-09-13</td><td>3.7</td><td>2.0 ± 0.8</td></tr><tr><td>NICER</td><td>3033290126</td><td>2020-09-15</td><td>1.9</td><td>1.5 ± 0.4</td></tr><tr><td>NuSTAR</td><td>80602315006</td><td>2020-09-17</td><td>78.1</td><td>0.8 ± 0.1</td></tr><tr><td>XMM-Newton</td><td>0872990101</td><td>2020-09-17</td><td>9.9</td><td>1.4 ± 0.2</td></tr></tbody></table>

<table><caption>**Table 5**
(Continued)</caption><thead><tr><th>Telescope</th><th>Obs ID</th><th>Start Day (YYYY-MM-DD)</th><th>Exposure<sup>a</sup> (ks)</th><th>Pulsed Flux<sup>b</sup> (10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>NICER</td><td>3033290127</td><td>2020-09-19</td><td>2.3</td><td>1.1 ± 0.3</td></tr><tr><td>NICER</td><td>3033290128</td><td>2020-09-22</td><td>2.4</td><td>1.4 ± 0.3</td></tr><tr><td>NICER</td><td>3033290130</td><td>2020-09-30</td><td>2.4</td><td>1.0 ± 0.4</td></tr><tr><td>XMM-Newton</td><td>0872990301</td><td>2020-10-02</td><td>19.1</td><td>1.1 ± 0.1</td></tr><tr><td>NICER</td><td>3033290131</td><td>2020-10-03</td><td>2.9</td><td>1.2 ± 0.3</td></tr><tr><td>NICER</td><td>3033290132</td><td>2020-10-05</td><td>1.8</td><td>1.2 ± 0.4</td></tr><tr><td>NICER</td><td>3033290133</td><td>2020-10-07</td><td>1.0</td><td>2.0 ± 0.6</td></tr><tr><td>NICER</td><td>3033290134</td><td>2020-10-07</td><td>0.6</td><td>1.8 ± 0.6</td></tr><tr><td>NuSTAR</td><td>80602315008</td><td>2020-10-09</td><td>153.4</td><td>0.7 ± 0.1</td></tr><tr><td>NICER</td><td>3033290135</td><td>2020-10-11</td><td>2.8</td><td>1.4 ± 0.3</td></tr><tr><td>NICER</td><td>3033290136</td><td>2020-10-13</td><td>1.8</td><td>1.3 ± 0.4</td></tr><tr><td>NICER</td><td>3033290137</td><td>2020-10-19</td><td>3.4</td><td>1.1 ± 0.3</td></tr><tr><td>XMM-Newton</td><td>0872990401</td><td>2020-10-23</td><td>33.9</td><td>1.1 ± 0.08</td></tr><tr><td>NICER</td><td>3033290138</td><td>2020-10-27</td><td>1.5</td><td>1.0 ± 0.5</td></tr><tr><td>NICER</td><td>3033290139</td><td>2020-10-27</td><td>1.5</td><td>1.1 ± 0.5</td></tr><tr><td>NICER</td><td>3033290140</td><td>2020-10-29</td><td>2.6</td><td>1.0 ± 0.4</td></tr><tr><td>NICER</td><td>3033290141</td><td>2020-10-31</td><td>2.2</td><td>1.1 ± 0.4</td></tr><tr><td>NICER</td><td>3033290142</td><td>2020-11-02</td><td>1.7</td><td>0.7 ± 0.4</td></tr><tr><td>NICER</td><td>3033290143</td><td>2020-11-04</td><td>3.1</td><td>0.7 ± 0.3</td></tr><tr><td>NICER</td><td>3033290144</td><td>2020-11-04</td><td>3.1</td><td>1.1 ± 0.4</td></tr><tr><td>NICER</td><td>3033290146</td><td>2020-11-09</td><td>3.0</td><td>0.9 ± 0.3</td></tr><tr><td>NICER</td><td>3033290147</td><td>2020-11-11</td><td>3.6</td><td>&lt;0.5</td></tr><tr><td>NICER</td><td>3033290148</td><td>2020-11-13</td><td>1.9</td><td>0.9 ± 0.4</td></tr><tr><td>NICER</td><td>3033290149</td><td>2020-11-14</td><td>1.9</td><td>&lt;1</td></tr><tr><td>NICER</td><td>3598011201</td><td>2020-11-18</td><td>1.0</td><td>0.9 ± 0.2</td></tr><tr><td>NICER</td><td>3598011202</td><td>2020-11-21</td><td>3.1</td><td>0.7 ± 0.3</td></tr><tr><td>NICER</td><td>4607020101</td><td>2021-03-08</td><td>2.0</td><td>1.1 ± 0.3</td></tr><tr><td>NICER</td><td>4607020201</td><td>2021-03-09</td><td>2.4</td><td>0.9 ± 0.4</td></tr><tr><td>NICER</td><td>4607020203</td><td>2021-04-14</td><td>7.9</td><td>0.6 ± 0.2</td></tr><tr><td>NICER</td><td>4607020204</td><td>2021-04-15</td><td>4.2</td><td>0.9 ± 0.3</td></tr><tr><td>NICER</td><td>4607020301</td><td>2021-03-11</td><td>3.0</td><td>0.9 ± 0.3</td></tr><tr><td>NICER</td><td>4607020401</td><td>2021-04-17</td><td>2.0</td><td>0.6 ± 0.4</td></tr><tr><td>NICER</td><td>4607020501</td><td>2021-05-03</td><td>7.9</td><td>0.6 ± 0.2</td></tr><tr><td>NICER</td><td>4607020601</td><td>2021-06-08</td><td>7.5</td><td>0.6 ± 0.2</td></tr><tr><td>NICER</td><td>4607020701</td><td>2021-07-14</td><td>7.0</td><td>0.5 ± 0.2</td></tr><tr><td>NICER</td><td>4607020802</td><td>2021-07-23</td><td>6.6</td><td>0.5 ± 0.2</td></tr><tr><td>NICER</td><td>4607020804</td><td>2021-08-09</td><td>3.6</td><td>0.5 ± 0.3</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Including only the Good Time Interval.

<sup>b</sup> Excluding PWN and SNR emission and in the energy range of 1–12 keV.

(This table is available in machine-readable form in the online article.)

Our spectral analysis indicates that there exist at least two distinct components contributing to the pulsed spectrum (over the 1.0–70 keV energy range): a single-temperature blackbody and a power law, both of which are absorbed by intervening hydrogen gas along the line of sight. Careful spectral modeling implies the gradual fading of the thermal component during the outburst, while the power law remains largely unchanged, which was also noted recently by C.-P. Hu et al. (2023). We illustrate the evolution in the best-fit blackbody temperatures and radii during the outburst of the source in Figure 10, comparing these values with that derived in C.-P. Hu et al. (2023), finding good agreement between the two studies in general. The temperature of the thermal component seems to correlate with the blackbody (and total) flux, showing an enhancement around MJD 59,082, while gradually decaying to
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The evolution in the best-fit blackbody temperatures and radii (top panels), and the pulsed flux and total flux (bottom panels) during the 2020 outburst of PSR J1846−0258. The values derived in this study (black points) are compared with those derived in C.-P. Hu et al. (2023, blue points). The upper limits (quoted at the $1\sigma$ confidence level) are indicated by the arrows.</figcaption></figure>

pre-outburst values a couple of months later. The derived blackbody radius appears constant for large durations of the outburst between 0.8 and 1.2 km, before decreasing to ~0.3–0.5 km at late stages of the outburst (i.e., after MJD 59,122).

The gradual fading of the thermal component could qualitatively be explained by the model of A. M. Beloborodov (2009), due to the dissipation of the twisted bundle of field lines in the magnetosphere, such that a smaller number of charged particles impact upon the NS surface. When the source returns to quiescence, the power law begins to dominate the pulsed emission (see, e.g., L. Kuiper et al. 2018 and E. V. Gotthelf et al. 2021 for an analysis over a broader energy range).

Moreover, we estimate that the total energy released during the outburst (assuming an exponential model for the flux decay, isotropic emission in the 1–12 keV energy band, and a source distance of 5.8 kpc) is ~1.4 × $10^{42}$ ergs. This is within a factor of 3 of the same quantity found for the previous outburst (in the 2–10 keV band) by F. P. Gavrilov et al. (2008). In comparison with other magnetar outbursts (see Figure 6 in F. Coti Zelati et al. 2018) PSR J1846−0258 does not show any dissimilarities with canonical magnetars despite its hybrid nature.

## Acknowledgments

R.S., N.R., and A.I. are supported by the H2020 ERC Consolidator Grant "MAGNESIA" under grant agreement No. 817661. F.C.Z. is supported by a Ramón y Cajal fellowship (grant agreement RYC2021-030888-I). R.S. also acknowledges financial support from the Italian Ministry for University and Research, through the grants 2022Y2T94C (SEAWIND) and from INAF through LG 2023 BLOSSOM. This work was also partially supported by the program Unidad de Excelencia María

de Maeztu CEX2020-001058-M. We thank Emilie Parent for useful discussions on the timing analysis. We also acknowledge partial support from grant SGR2021-01269 (PI: Graber/Rea). A.I.'s work has been carried out within the framework of the doctoral program in Physics at the Universitat Autonoma de Barcelona. Part of the research activities described in this paper were carried out with contribution of the NextGenerationEU funds within the National Recovery and Resilience Plan (PNRR), Mission 4—Education and Research, Component 2—From Research to Business (M4C2), Investment Line 3.1—Strengthening and creation of Research Infrastructures, Project IR0000026-Next Generation Croce del Nord. The Parkes radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42), which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri people as the Traditional Owners of the Observatory site. This work has made use of the ` matplotlib ` (J. D. Hunter 2007), ` numpy ` (C. R. Harris et al. 2020), and ` scipy ` (P. Virtanen et al. 2020) libraries in ` pythonv3.8`.

## ORCID iDs

R. Sathyaprakash https://orcid.org/0000-0002-5254-3969

N. Rea https://orcid.org/0000-0003-2177-6388

F. Coti Zelati https://orcid.org/0000-0001-7611-1581

A. Borghese https://orcid.org/0000-0001-8785-5922

M. Pilia https://orcid.org/0000-0001-7397-8091

M. Trudu https://orcid.org/0000-0002-1530-0474

M. Burgay https://orcid.org/0000-0002-8265-4344

R. Turolla https://orcid.org/0000-0003-3977-8760

S. Zane https://orcid.org/0000-0001-5326-880X
P. Esposito https://orcid.org/0000-0003-4849-5092

S. Mereghetti https://orcid.org/0000-0003-3259-7801

S. Campana https://orcid.org/0000-0001-6278-1576

D. Götz https://orcid.org/0000-0001-9494-0981

A. Y. Ibrahim https://orcid.org/0000-0002-5663-1712

G. L. Israel https://orcid.org/0000-0001-5480-6438

A. Possenti https://orcid.org/0000-0001-5902-3731

A. Tiengo https://orcid.org/0000-0002-6038-1090

## References

Alpar, M. A., Cheng, K. S., & Pines, D. 1989, ApJ, 346, 823

Archibald, R. F., Kaspi, V. M., Beardmore, A. P., Gehrels, N., & Kennea, J. A. 2015, ApJ, 810, 67

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

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

Ashton, G., Jones, D. I., & Prix, R. 2017, MNRAS, 467, 164

Beloborodov, A. M. 2009, ApJ, 703, 1044

Blumer, H., Safi-Harb, S., McLaughlin, M. A., & Fiore, W. 2021, ApJL, 911, L6

Brook, P. R., Karastergiou, A., Johnston, S., et al. 2016, MNRAS, 456, 1374

Camilo, F., Ransom, S. M., Halpern, J. P., et al. 2006, Natur, 442, 892

Cordes, J. M., & McLaughlin, M. A. 2003, ApJ, 596, 1142

Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961

Esposito, P., Rea, N., & Israel, G. L. 2021, in Magnetars: A Short Review and Some Sparse Considerations, ed. T. M. Belloni, M. Méndez, & C. Zhang (Berlin: Springer), 97

Freire, P. C. C., & Ridolfi, A. 2018, MNRAS, 476, 4794

Gavriil, F. P., Gonzalez, M. E., Gotthelf, E. V., et al. 2008, Sci, 319, 1802

Gotthelf, E. V., Safi-Harb, S., Straal, S. M., & Gelfand, J. D. 2021, ApJ, 908, 212

Gotthelf, E. V., Vasisht, G., Boylan-Kolchin, M., & Torii, K. 2000, ApJL, 542, L37

Güercinoğlu, E., Köksal, E., & Güver, T. 2023, MNRAS, 518, 5734

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

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

Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1

Helfand, D. J., Collins, B. F., & Gotthelf, E. V. 2003, ApJ, 582, 783

Hobbs, G. 2009, arXiv:0911.5534

Hobbs, G., Lyne, A. G., & Kramer, M. 2010, MNRAS, 402, 1027

Hobbs, G., Manchester, R. N., Dunning, A., et al. 2020, PASA, 37, e012

Hu, C.-P., Kuiper, L., Harding, A. K., et al. 2023, ApJ, 952, 120

Hunter, J. D. 2007, CSE, 9, 90

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

Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151

Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261

Krimm, H. A., Lien, A. Y., Page, K. L., Palmer, D. M., & Tohuvavohu, A. 2020, ATel, 13913, 1

Kuiper, L., Harding, A. K., Enoto, T., et al. 2020, ATel, 13985, 1

Kuiper, L., & Hermsen, W. 2009a, A&A, 501, 1031

Kuiper, L., & Hermsen, W. 2009b, ATel, 2151, 1

Kuiper, L., Hermsen, W., & Dekker, A. 2018, MNRAS, 475, 1238

Kumar, H. S., & Safi-Harb, S. 2008, ApJL, 678, L43

Kumar, P., Shannon, R. M., Flynn, C., et al. 2021, MNRAS, 500, 2525

Laha, S., Barthelmy, S. D., Cummings, J. R., et al. 2020, GCN, 28193, 1

Leahy, D. A., & Tian, W. W. 2008, A&A, 480, L25

Livingstone, M. A., Kaspi, V. M., Gotthelf, E. V., & Kuiper, L. 2006, ApJ, 647, 1286

Livingstone, M. A., Ng, C. Y., Kaspi, V. M., Gavriil, F. P., & Gotthelf, E. V. 2011, ApJ, 730, 66

Lorimer, D., & Kramer, M. 2004, Handbook of Pulsar Astronomy, Cambridge Observing Handbooks for Research Astronomers, Vol. 4 (Cambridge: Cambridge Univ. Press)

Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

Majid, W. A., Pearlman, A. B., Prince, T. A., et al. 2020, ATel, 13988, 1

McLaughlin, M. A., Lyne, A. G., Keane, E. F., et al. 2009, MNRAS, 400, 1431

Mickaliger, M. B., Rajwade, K., Stappers, B., Lyne, A., & Preston, L. L. 2020, ATel, 13950, 1

NASA High Energy Astrophysics Science Archive Research Center (HEASARC) 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Nice, D., Demorest, P., Stairs, I., et al. 2015, Tempo: Pulsar Timing Data Analysis, Astrophysics Source Code Library, ascl:1509.002

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, JMLR, 12, 2825

Rajwade, K. M., Stappers, B. W., Lyne, A. G., et al. 2022, MNRAS, 512, 1687

Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning (Cambridge, MA: MIT Press)

Rea, N., Esposito, P., Turolla, R. et al. 2010, Sci, 330, 944

Straal, S. M., Gelfand, J. D., & Eagle, J. L. 2023, ApJ, 942, 103

Terlouw, J. P., & Vogelaar, M. G. R. 2014, Kapteyn Package, version 3.4, Astronomical Source Code Library, ascl:1611.010

Uzuner, M., Keskin, Ö., Kaneko, Y., et al. 2023, ApJ, 942, 8

Verbiest, J. P. W., Weisberg, J. M., Chael, A. A., Lee, K. J., & Lorimer, D. R. 2012, ApJ, 755, 39

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wang, F. Y., Zhang, G. Q., Dai, Z. G., & Cheng, K. S. 2022, NatCo, 13, 4382

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914