Long-term Coherent Timing of the Accreting Millisecond Pulsar IGR J17062–6143

Between 2017 August and 2020 August, NICER observed IGR J17062 for a total unfiltered exposure of 372 ks. Based on their spacing in time, these data are naturally grouped into 10 distinct epochs. We list these groups and their respective ObsIDs in Table 1.

We processed all data using nicerdas version 7a, which was released with heasoft version 6.27.2. By default, this pipeline selects only those epochs with a pointing offset <54'', a bright-Earth limb angle >30°, a dark-Earth limb angle >15°, and which are outside of the South Atlantic Anomaly. It additionally screens for signs of increased background emission by retaining only those data epochs that were collected at times when the rate of reset triggers (undershoots) is <200 ct/s/detector and the rate of high-energy events (overshoots) is <1 ct/s/detector and <1.52 × COR_SAX^−0.633, where cor_sax gives the cutoff rigidity of Earth's magnetic field. We find that these latter two criteria are often too conservative when applied to the observations of IGR J17062. Following Bult et al. (2020a), we first smoothed the overshoot light curve using 5 s windows to reduce noise and then only retained those epochs when the absolute rate was less than 1.5 ct/s/detector and increased the scaling factor of the cor_sax expression filter from 1.52 to 2.0. The undershoot filter was only an issue in data groups 3, 4, 5, and 8. These data were all collected at times when the Sun angle was comparatively low (<65°), which increases the optical load on the instrument and affects the low-energy background. We increased the undershoot rate filter for these data from 200 ct/s/detector to 400 ct/s/detector. For all ObsIDs, we compare the light curves obtained with default screening criteria to those obtained with our relaxed parameters to ensure no spurious features are introduced into the data selection.

Inspecting the light curves, we find that a number of individual exposures in the third data group show near identical light curve profiles in which the overall count-rate drops linearly from about 40 ct s ⁻¹ to 0 ct s⁻¹. These decays coincide with a decrease in the number of stars traced by the star tracker and reflect an obstruction of the field of view rather than being of astrophysical origin. ¹⁰ We reprocessed these particular ObsID with the additional requirement that the number of stars in the star tracker (st_stars) be larger than 38, which effectively removed all spurious trends.

After screening our data, we are left with 202 ks of clean exposure (as opposed to the 187 ks retained under standard screening criteria). We applied barycenter corrections to these cleaned data using the Swift/UVOT position of Ricci et al. (2008) and the JPL-DE405 solar system ephemeris (Standish 1998). All dates reported in this work are therefore in terms of Barycentric Dynamic Time (TDB).

XMM-Newton observed IGR J17062 on 2016 September 13–15 (MJD 57645; see van den Eijnden et al. 2018 for a detailed analysis). These data were collected with the epic-pn camera in timing mode, yielding a relative time resolution of 30 μs. The absolute timing calibration of XMM-Newton has an accuracy of 100 μs (Rosen 2020), which translates to a systematic phase uncertainty of 0.02 cycles for IGR J17062. Hence, these data can be used to extend the coherent timing analysis of our NICER observations.

We processed the XMM-Newton data in SAS v18 using the latest version of the calibration files. Standard screening filters were applied, i.e, we kept only those events with photon energies in the 0.4–10 keV range, with PATTERN ≤ 4 and screening FLAG = 0. We extracted the source event list from rawx columns [34:42] and obtained background events from rawx columns [6:14]. The source light curve showed a mean count rate of 18 ct s⁻¹ over the 61 ks exposure, with a background contribution of 0.26 ct s⁻¹. Finally, we used barycen to apply barycentric corrections to our source data, based on the same source coordinates and ephemerides used for NICER.

RXTE observed IGR J17062 on 2008 May 3 (MJD 54589.0) for a total good time exposure of 1.1 ks in the proportional counter array (PCA). The data were collected with the PCA operating in event mode, using 64 energy channels and a time resolution of 1/8192 s. A detailed analysis of these data is presented by Strohmayer & Keek (2017). Here, we applied standard processing and screening methods to extract a photon event list for this observation and used faxbary to apply barycentering corrections. The background rate was estimated using pcabackest.

Following the approach of Strohmayer et al. (2018), we separately search each data group for the presence of 163 Hz pulsations by optimizing the ${Z}_{1}^{2}$ score (Buccheri et al. 1983) as a function of spin frequency and the binary orbit's time of ascending node, T_asc. We adopt the ephemeris reported by Strohmayer et al. (2018) as our trial solution ¹² and correct the photon arrival times to remove the Doppler modulation imposed by the circular binary orbit. We then compute the ${Z}_{1}^{2}$ score as

$$\begin{eqnarray}&&{Z}_{1}^{2}=\displaystyle \frac{2}{N}\left[{\left(\displaystyle \sum _{j=1}^{N}\cos {\varphi }_{j}\right)}^{2}+{\left(\displaystyle \sum _{j=1}^{N}\sin {\varphi }_{j}\right)}^{2}\right],\end{eqnarray} \tag{ 1 }$$

with φ_j = 2π ν_i t_j , where the t_j give the list of N photon arrival times in a group, and ν_i is the test frequency. We evaluate this score on a grid of trial frequencies spanning ±40 μHz around the source spin frequency. The window width was conservatively chosen such that we are sensitive to a long-term spin-frequency derivative smaller than 4 × 10⁻¹³ Hz s⁻¹ over the three-year span between the first and last data groups. The grid resolution is set to 1/(8T), where T is the duration of the respective data group, that is, the time between the first and last photon events in that group. For T_asc we adopt a grid spanning one orbital period that we center on the middle of the respective data group's time interval. This grid has a resolution of 25 s, which is equivalent to 4° orbital longitude.

We applied the ${Z}_{1}^{2}$ optimization method to all data groups, excluding group 9 which is dominated by the X-ray burst response, and groups 4 and 5, which do not have sufficient exposure for this type of broadly defined search. Coherent 163 Hz pulsations were recovered in all searched groups. We determined uncertainties on the resulting spin frequency and T_asc measurements by scanning the ${Z}_{1}^{2}$ space for the boundary interval where ${\rm{\Delta }}{Z}_{1}^{2}=9$, which corresponds to the 3σ confidence boundaries (Markwardt et al. 2002). To further verify the validity of these parameter uncertainties, we simulated 1000 Poisson sampled event lists from the timing model while maintaining the same count rates and good time intervals as our observations. Running the ${Z}_{1}^{2}$ optimization method on these simulated data, we obtained parameter uncertainties that were entirely consistent with those found through the ${\rm{\Delta }}{Z}_{1}^{2}=9$ interval scan. The resulting measurements and uncertainties are shown in Figure 2, with the detailed parameter values reported in the appendix.

Considering the pulse-frequency measurements, we observe a modest linear drift across the data groups (Figure 2, top panel). To verify if this drift is statistically significant, we first fit the measurements using a constant frequency model. We obtain a best-fit χ² of 39 for six degrees of freedom, yielding a p value of 7 × 10⁻⁷. Hence, this fit firmly rejects the constant spin-frequency hypothesis. Fitting the measurements with a linear function instead, we obtain an acceptable fit, if slightly underdispersed at a χ² of 1.0 for five degrees of freedom. The best-fit parameters of this linear fit give a spin-frequency measurement of ν = 163.65611072 ± 4 × 10⁻⁸ Hz at a reference time of 58522.3 TDB, with a significantly detected long-term spin-up of $\dot{\nu }=(5.3\pm 0.9)\times {10}^{-15}\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

Considering the T_asc measurements, we observe a similar long-term evolution in our measurements. In Figure 2 (bottom panel), we plot the measured T_asc relative to the predicted values based on the constant orbital period model of Strohmayer et al. (2018). These measurements are well described using a linear function, yielding a best-fit χ² of 8.4 for five degrees of freedom. The resulting best-fit orbital period is P_b = 2278.2110 s ± 0.3 ms, which is 3 ms larger than the orbital period reported in Strohmayer et al. (2018) and well within their 12 ms uncertainty region. Adding a quadratic term to the model, we can further place a 95% confidence upper limit on the long-term derivative of the orbital period of $| {\dot{P}}_{b}| \lt {10}^{-11}\,{\rm{s}}\,{{\rm{s}}}^{-1}$.

Using the spin-frequency and binary-period parameters obtained from these linear models, we refined the timing solution and extrapolated the ephemeris to the epoch of the XMM-Newton observation. We then selected the 0.4–6 keV events and corrected their arrival times to the pulsar's reference frame using this predicted ephemeris. Calculating the ${Z}_{1}^{2}$ score for this data, we obtain a score of ${Z}_{1}^{2}=18,$ which has an associated single-trial p value of 1 × 10⁻⁴. Hence, the pulsations are significantly detected in the XMM-Newton observation at the 4σ level.

We also used the interim timing solution to consider NICER data groups 4, 5, and 9. In each case, we corrected the arrival times for the predicted binary orbit and determined the ${Z}_{1}^{2}$ score. Significant pulsations we recovered in groups 4 (p value of 3 × 10⁻⁴; 3.6σ) and 5 (p value of 7 × 10⁻⁶; 4.5σ), but not in group 9 (p value of 0.77). Because the source luminosity varies drastically across data group 9, we also searched each individual ObsID; however, this did not yield a pulse detection either. Finally, we attempted to suppress the influence of the X-ray burst by processing the data in reverse order. That is, we first searched the last ObsID for pulsations. As none were detected, we added the second to last and-so forth. At no point in this process did the pulse significance exceed 1σ. In the absence of a detection, we calculate a 95% confidence level upper limit on the pulse fraction of 0.8% if we include all ObsIDs in this group. We note, though, that this limit relaxes to 1.1% if exclude the first four ObsIDs, which show the largest swing in source intensity (see Figure 1).

A fully coherent description of the pulse signal requires that our timing solution is precise enough to exactly predict the number of orbital revolutions and stellar rotations that take place during the time intervals between the data groups. Adopting the interim timing solution obtained in the previous section, we find that the orbital phase is readily extrapolated across the time span of our data (as is clear from Figure 2). The stellar rotation, on the other hand, requires more attention. The uncertainty in the predicted pulse phase grows approximately as

$$\begin{eqnarray}&&{\sigma }_{\phi }^{2}(t)={\left({\sigma }_{\nu }t\right)}^{2}+{\left(\displaystyle \frac{1}{2}{\sigma }_{\dot{\nu }}{t}^{2}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where σ_ν gives the uncertainty on the spin frequency and ${\sigma }_{\dot{\nu }}$ gives the uncertainty on its derivative. Using the uncertainties obtained in Section 4.1, we find that the timing solution loses coherence (σ_ϕ exceeds 0.5) when t ≈ 200 days. Hence, we can coherently connect some data groups, but not all of them at once. To further refine the timing solution, we therefore start by combining data groups 2, 3, and 4, which have the shortest time intervals between them. For each group, we fold the data to a pulse profile, which we fit with a sinusoid to measure the pulse-phase residual relative to the model. The set of phase residuals are then modeled using a quadratic function, so that we obtain refined measurements for ν and $\dot{\nu }$, thus increasing the coherence time. By iteratively adding the nearest data group to the analysis, we gradually converge to a single timing solution that coherently describes all data. This approach yielded a final spin frequency of ν = 163.6561106613 ± 1.0 × 10⁻⁹ Hz and a derivative of $\dot{\nu }=(3.73\pm 0.07)\times {10}^{-15}\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$ measured relative to a spin epoch of 58522.3 TDB.

While the iterative approach allowed us to determine a coherent solution for the whole data set, this method is at risk of converging to a local minimum if the timing solution is not unique. To verify the uniqueness of the solution, we took a numerical approach. We scanned through the joint 5σ confidence region in ν and $\dot{\nu }$ while taking steps of 10⁻¹⁰ Hz in frequency space and steps of 10⁻¹⁸ Hz s⁻¹ in the derivative. For each trial, we compute the phase residuals for the XMM-Newton observation and each of the NICER data groups (excluding only group 9) and calculate the χ² of this model. As shown in Figure 3, this χ² scan yielded only one plausible timing solution, which had the same $(\nu ,\dot{\nu })$ parameters found previously through the iterative approach.

Zoom In Zoom Out Reset image size

We adopted the spin frequency and spin-frequency derivative measurements associated with the minimum in χ² space and constructed the phase residuals. We then fit these residuals jointly for the frequency evolution and the binary-orbit parameters (orbital period, projected semimajor axis, and T_asc). This fit yielded a slightly improved description of the pulse arrival times and higher-precision measurements of the orbital parameters. We further attempted adding a second frequency derivative to the model, but this did not improve the fit. Including orbital eccentricity or a binary-period derivative was not required either. Finally, optimizing the pulse-phase residuals as a function of the source coordinate did not improve the residuals either, suggesting the source position is correct to within a 90% confidence uncertainty region of 0.1''. The complete timing model is listed in Table 2, while the resulting phase residuals are shown in Figure 4.

Zoom In Zoom Out Reset image size

Table 2. NICER Timing Solution

Parameter	Value	Uncertainty
R.A. (J2000)	17:06:16.29	01
Decl. (J2000)	–61:42:40.6	01
ν (Hz)	163.6561106623	1.0 × 10−9
$\dot{\nu }$ (Hz s −1)	3.77 × 10−15	9 × 10−17
$| \ddot{\nu }|$ (Hz s −2)	<1.3 × 10−24
Spin epoch	58522.3
Pb (s)	2278.21124	2 × 10−5
${\dot{P}}_{b}$ (s s −1)	8.4 × 10−12	2 × 10−12
Semimajor axis (lt-s)	0.003963	6 × 10−6
Tasc (TDB)	58588.78464	5 × 10−5
Eccentricity	<0.03

Note. Source coordinates were adopted from Ricci et al. (2008), with improved uncertainties obtained through the timing model. Upper limits are quoted at 95% confidence.

Download table as:  ASCIITypeset image

Using the refined timing solution, we can revisit the 2008 RXTE observation to test for consistency. Extrapolating the spin frequency and T_asc to the RXTE epoch (54589.5 TDB) and folding these data, we measure a pulse fraction of about 7%. Given that Strohmayer & Keek (2017) report a pulse fraction of 9%, this result indicates our timing solution is still suboptimal. We therefore optimized the local timing solution by again determining the peak ${Z}_{1}^{2}$ score as a function of the time of ascending node and the spin frequency, finding T_asc = 54589.5510 ± 0.0003 TDB and ν = 163.65612 ± 5 × 10⁻⁵ Hz. This frequency measurement is consistent with the long-term frequency solution, deviating only by 0.1σ. The time of ascending node, on the other hand, arrives 233 s later than predicted, which constitutes a 3.3σ shift. Combining this measurement with the T_asc values obtained from NICER and XMM-Newton, we find that the long-term drift in T_asc is incompatible with a linear trend (χ²/dof = 25.9/7) but well described by a quadratic function (χ²/dof = 8.9/6). The binary-period derivative obtained through this fit is ${\dot{P}}_{b}=(8.4\pm 2.0)\times {10}^{-12}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, indicating that the binary orbit is expanding. The reference T_asc and orbital period obtained through this fit are consistent with those found previously (to within 1σ).

Considering the pulse amplitudes across the different data groups, we observe a clear evolution in pulse fraction over time. This trend could be well described as a sinusoid (χ²/dof = 6.4/6), yielding an amplitude of (1.25 ± 0.16)% over a mean pulse fraction of (1.94 ± 0.08)%, with a period measurement of 1210 ± 40 days, or about 3.3 years (see Figure 4). The 1–10 keV X-ray flux of these data groups does not exhibit any such oscillation, instead showing a random 13% rms scatter around a mean flux of ≈6 × 10⁻¹¹ erg s⁻¹ cm⁻².

To investigate the energy dependence of the pulsations, we folded all NICER data to pulse profiles in several energy bins spanning 0.4–10 keV. The bins were constructed adaptively to be multiples of 0.5 keV, such that each bin contained at least 10⁴ counts and the pulsation was detected at >5σ significance. As shown in Figure 5, we find that the time-averaged fractional pulse amplitude increases with energy to a peak of about 4.5% between 5 and 10 keV. For the pulse phase, we find that the pulsations have a roughly constant phase below 2 keV, but show a decreasing phase lag toward higher energies, i.e., the softer photons (<2 keV) tend to arrive after the harder photons (>2 keV).

Zoom In Zoom Out Reset image size

For comparison, we also calculated the energy dependence of the pulse detection in the 2008 RXTE observation. Using the local timing solution, we adaptively binned individual event channels such that each bin contained a pulse detection of at least >5σ. We find the pulsations could be resolved in five separate energy bins, which are shown in Figure 5 using red diamonds.

We have measured the long-term spin-frequency derivative of IGR J17062 to be $\dot{\nu }=(3.73\pm 0.03)\times {10}^{-15}\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$, indicating the neutron star is slowly spinning up under the influence of sustained mass accretion onto the star. Positive pulse-frequency derivatives have been previously reported for several other AMXPs during outburst (Di Salvo & Sanna 2020), with magnitudes on the order of 10⁻¹⁴ ∼ 10⁻¹³ Hz s⁻¹. In all cases, these derivatives were measured by modeling the quadratic drift in pulse phase. However, flux-dependent timing noise can bias such phase measurements, making it difficult to disentangle what part of the observed phase drift is intrinsic to the pulsar rotation and what part is instead due to accretion-driven variability in the hot spot (Patruno et al. 2009; Bult et al. 2020a). Our spin-frequency derivative measurement for IGR J17062 does not suffer from this ambiguity: the spin change is directly visible in frequency space.

In the conventional picture of accretion onto a magnetized star, the spin evolution of the star is governed by the interaction between the stellar magnetosphere and the inner accretion disk. The stellar magnetic field truncates the accretion disk at r_m, the magnetospheric radius. Just outside r_m, the field threads the disk, leading to an exchange of angular momentum between the star and the accretion flow that depends on the difference in their angular velocities. While the precise torques acting on the star depend on the specifics of the interaction model, the torque is generally expressed as

$$\begin{eqnarray}&&{N}_{\mathrm{total}}={\dot{M}}_{\mathrm{ns}}\sqrt{{{GM}}_{\mathrm{ns}}{r}_{{\rm{m}}}}+{N}_{\mathrm{field}},\end{eqnarray} \tag{ 3 }$$

where ${\dot{M}}_{\mathrm{ns}}$ is the mass accretion rate onto the neutron star, G the gravitational constant, and M_ns the mass of the neutron star. The N_field term expresses the torque exchange between the magnetosphere and the inner accretion disk, and depends on specific assumptions about the microphysics of the magnetic coupling. Rappaport et al. (2004) studied the specific case of an accretion disk around a millisecond pulsar and derived

$$\begin{eqnarray*}{N}_{\mathrm{field}}=\left\{\begin{array}{ll}\tfrac{{\mu }^{2}}{3{r}_{{\rm{c}}}^{3}}\left[\tfrac{2}{3}-2{\left(\tfrac{{r}_{{\rm{c}}}}{{r}_{{\rm{m}}}}\right)}^{3/2}+{\left(\tfrac{{r}_{{\rm{c}}}}{{r}_{{\rm{m}}}}\right)}^{3}\right] & {r}_{{\rm{m}}}\lt {r}_{{\rm{c}}},\\ -\tfrac{{\mu }^{2}}{9{r}_{{\rm{c}}}^{3}} & {r}_{{\rm{m}}}\gt {r}_{{\rm{c}}}\end{array}\right.,\end{eqnarray*}$$

where μ is the magnetic dipole moment and

$$\begin{eqnarray}&&{r}_{{\rm{c}}}={\left(\displaystyle \frac{{{GM}}_{\mathrm{ns}}}{{{\rm{\Omega }}}_{s}^{2}}\right)}^{1/3},\end{eqnarray} \tag{ 4 }$$

the corotation radius, with Ω_s = 2π ν the angular velocity of the neutron star. Adopting a 1.4M_⊙ neutron star, we find that the corotation radius of IGR J17062 is 56 km.

The measured spin-frequency derivative gives us a measure of the total torque acting on the star as ${N}_{\mathrm{total}}=2\pi I\dot{\nu }$, where I = (0.5–2) × 10⁴⁵ g cm² gives the neutron star moment of inertia (Friedman et al. 1986; Steiner et al. 2015). Substituting this expression into Equation (3), we can place a lower limit on the magnetic dipole moment of μ ≥ 1 × 10²⁶ G cm³, with a magnetospheric radius of r_m ≥ 50 km. We can obtain more precise measurements if we further assume that the magnetospheric radius scales with the mass accretion rate as

$$\begin{eqnarray}&&{r}_{{\rm{m}}}=\xi {\left(\displaystyle \frac{{\mu }^{4}}{{{GM}}_{\mathrm{ns}}{\dot{M}}^{2}}\right)}^{1/7},\end{eqnarray} \tag{ 5 }$$

with scaling factor ξ ≈ 0.5 (Long et al. 2005; Bessolaz et al. 2008; Zanni & Ferreira 2013). Assuming the X-ray flux is a good tracer of the mass accretion rate, we adopt the 0.3–79 keV luminosity measured by van den Eijnden et al. (2018) to find that the accretion rate onto the neutron star is ${\dot{M}}_{\mathrm{ns}}=2.5\times {10}^{-11}{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Considering the energetics and recurrence time of the X-ray bursts, Keek et al. (2017) derive a similar accretion rate, suggesting that this estimate is quite robust. Solving the set of equations gives μ = 2.9 × 10²⁶ G cm³ with a truncation radius of 52 km.

Obviously, the derived dipole moment is highly model dependent, as somewhat different prescriptions of the disk–field interaction will produce slightly different outcomes. There may also be additional torques acting on the neutron star, such as the torque associated with a propeller mechanism (Illarionov & Sunyaev 1975). There is reason to suspect that a propeller mechanism might be active in IGR J17062; for one, the mass flow rate through the accretion disk has been estimated to be 1.8 × 10⁻¹⁰ M_⊙ yr⁻¹ (Hernández Santisteban et al. 2019), about an order of magnitude higher than the accretion rate onto the neutron star. Additionally, high-resolution X-ray spectroscopy (Degenaar et al. 2017b; van den Eijnden et al. 2018) revealed narrow emission lines, which can be interpreted as the blueshifted emission of an outflow. If we assume that IGR J17062 is in a weak propeller regime (Romanova et al. 2005; Ustyugova et al. 2006), then part of the accretion flow is ejected by the propeller, and the remainder accretes onto the neutron star. In this case, the material pressure exerted on the magnetosphere can be much greater than inferred from the X-ray flux, while the ejection mechanism itself applies a negative torque to the neutron star. The torque equation becomes

$$\begin{eqnarray}&&{N}_{\mathrm{total}}=\dot{M}\sqrt{{{GM}}_{\mathrm{ns}}{r}_{{\rm{m}}}}+{N}_{\mathrm{field}}+{N}_{\mathrm{propeller}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{N}_{\mathrm{propeller}}={\dot{M}}_{\mathrm{ej}}\sqrt{{{GM}}_{\mathrm{ns}}{r}_{{\rm{m}}}},\end{eqnarray} \tag{ 7 }$$

with $-{\dot{M}}_{\mathrm{ej}}={\dot{M}}_{\mathrm{disk}}-{\dot{M}}_{\mathrm{ns}}\lt 0$ expressing the rate at which mass is being propelled out of the system. If we assume that the magnetospheric radius depends on ${\dot{M}}_{\mathrm{disk}}$, we can again solve the equation to find μ = 6.3 × 10²⁶ G cm³ at a radius of 46 km.

Whether we include a propeller torque or not, we therefore find that the observed spin-up of the neutron star is consistent with a stellar magnetic field of B ≈ 5 × 10⁸ G and an accretion disk that is truncated close to the corotation radius. Here, it is interesting to note that the X-ray reflection spectroscopy of the Fe Kα emission line is subject to a degeneracy between the inner truncation radius of the disk and the system inclination (van den Eijnden et al. 2018). Hence, if we use the truncation radius obtained from the torque analysis to break this degeneracy, the system inclination would have to be 30°–35°.

By modeling the advancement in the time of ascending node passages between 2008 and 2020, we obtained a binary-period derivative measurement of ${\dot{P}}_{b}=(8.4\pm 2.0)\times {10}^{-12}\,{\rm{s}}\,{{\rm{s}}}^{-1}$. Hence, we find that the binary orbit is expanding with an evolutionary timescale of ${P}_{b}/{\dot{P}}_{b}=8.6\,\mathrm{Myr}$.

The apparent binary evolution of IGR J17062 is much faster than expected from theory. Generally, the mass transfer in an ultracompact binary is believed to be driven by angular momentum losses through gravitational wave emission (Kraft et al. 1962), causing the binary orbit to gradually widen with time. The expected binary-period derivative follows as (Rappaport et al. 1987; Verbunt 1993; di Salvo et al. 2008)

$$\begin{eqnarray}\begin{array}{rcl}{\dot{P}}_{b} & = & 4.4\times {10}^{-14}\left(\displaystyle \frac{{M}_{\mathrm{ns}}{M}_{c}}{{M}_{\odot }}\right){\left(\displaystyle \frac{{M}_{\mathrm{total}}}{{M}_{\odot }}\right)}^{-1/3}\\ & & \times \,{\left(\displaystyle \frac{{P}_{b}}{1\mathrm{hr}}\right)}^{-5/3}\displaystyle \frac{n-1/3}{n+5/3-2q}\,\,{\rm{s}}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 8 }$$

with M_c the mass of the companion star, q = M_c /M_ns the binary mass ratio, and n the index of the companion star's mass–radius relation (n = −1/3 for a degenerate star). Adopting a 1.4M_⊙ neutron star and a 30° inclination (see Section 6.1), we find a companion mass of M_c = 0.011M_⊙ from the mass function (Strohmayer et al. 2018) and obtain an expected orbital period derivative for IGR J17062 of 10⁻¹⁵ s s⁻¹, which is clearly inconsistent with our measurement.

Similarly, rapid orbital evolution has been observed in a number of other AMXPs and LMXBs (Patruno et al. 2017; Di Salvo & Sanna 2020). In a few particularly well-sampled sources, most notably SAX J1808.4–3658 (Bult et al. 2020a), the advancing T_asc shows complex residuals, suggesting that the orbital period evolution is modulated (quasi-)periodically over decades. In IGR J17062, such modulation is not apparent in the data. As the period-derivative measurement is strongly influenced by the data gap between the 2008 RXTE and 2016 XMM-Newton observations, however, we note that further monitoring is required to assess the stability of the period change.

Various models have been proposed to explain the anomalously rapid evolution observed in some X-ray binaries, most prominently donor mass loss, spin–orbit coupling, and enhanced magnetic braking (see, e.g., di Salvo et al. 2008; Burderi et al. 2009; Patruno et al. 2017; Sanna et al. 2017b, for in depth discussions). Each of these models, however, faces some challenges. Spin–orbit coupling relies on the presence of a variable mass quadrupole in the companion star (Applegate & Shaham 1994; Richman et al. 1994), however, a low-mass companion likely does not have a sufficient energy budget for this mechanism (Brinkworth et al. 2006), which is only exacerbated for IGR J17062, which has the smallest known mass function among stellar binaries (Strohmayer et al. 2018). Enhanced magnetic braking might play a role (Justham et al. 2006; Ginzburg & Quataert 2020), but only if the companion star can sustain a ∼1 kG magnetic field and is also losing mass through a stellar wind. In the following, we therefore limit our discussion to just the mass-loss model.

If we assume that mass transfer is nonconservative, then the matter being ejected from the binary system can carry away additional orbital angular momentum. Assuming that this mass ejection channel dominates the evolution timescale, the binary-period derivative is, to good approximation, given by the relation (Rappaport et al. 1987)

$$\begin{eqnarray*}&&\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}=-\displaystyle \frac{{\dot{M}}_{c}}{{M}_{c}},\end{eqnarray*}$$

so that the required mass-loss rate in the companion star must be ${\dot{M}}_{c}\approx 2\times {10}^{-9}{M}_{\odot }\,{\mathrm{yr}}^{-1}$.

As we observe an ongoing outburst from IGR J17062, at least some of the matter lost from the companion must be transferred to the neutron star. Following Rappaport et al. (1982, 1983), we can assume that a fraction (1 − β) of the matter lost by the companion is ejected from the system, carrying specific angular momentum α, expressed in units of the companion star's orbital angular momentum. The remainder of the mass loss is transferred to the neutron star/accretion disk, so we can write ${\dot{M}}_{\mathrm{disk}}=-\beta {\dot{M}}_{c}$, where the mass flow rate through the disk has been observationally constrained to ${\dot{M}}_{\mathrm{disk}}=1.8\times {10}^{-10}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Hernández Santisteban et al. 2019). Hence, about 90% of the mass lost by the donor will be immediately ejected, while the remaining 10% will pass through the disk and move toward the neutron star. Given that the mass accretion rate onto the neutron star is another order of magnitude lower (Keek et al. 2017; Strohmayer et al. 2018; van den Eijnden et al. 2018), it seems that most of the matter flowing through the disk will eventually also be ejected. This secondary outflow, however, likely occurs as a disk wind or a propeller, meaning that the specific orbital angular momentum it carries is that of the neutron star (${\alpha }_{\mathrm{ns}}\approx {({a}_{1}/{a}_{2})}^{2}=0.01$) rather than that of the companion star (α ≈ 1). Hence, the impact of the second outflow on the binary evolution is expected to be negligible. Under these assumptions, the mass-loss rate of the companion is given as (Rappaport et al. 1983; Verbunt 1993)

$$\begin{eqnarray}&&{\dot{M}}_{c}=\displaystyle \frac{32{G}^{3}}{5{c}^{5}}{\left(\displaystyle \frac{4{\pi }^{2}}{G}\right)}^{4/3}{\left(\displaystyle \frac{{M}_{\mathrm{ns}}}{{P}_{b}}\right)}^{8/3}\displaystyle \frac{{q}^{2}}{{\left(1+q\right)}^{1/3}{f}_{r}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{f}_{r}=\displaystyle \frac{5}{6}+\displaystyle \frac{n}{2}-\beta q-\displaystyle \frac{(1-\beta )(3\alpha +q)}{3(1+q)}.\end{eqnarray} \tag{ 10 }$$

Equating this expression to the mass-loss rate implied by the orbital period derivative, we can solve for α as a function of just the neutron star mass and system inclination (Figure 7). We find α = 0.72–0.91, which corresponds to an ejection radius whose specific orbital angular momentum lies somewhere between the specific orbital angular momentum carried by the outer edge of the accretion disk (∼0.75) and that of the inner Lagrange point (∼0.9).

Zoom In Zoom Out Reset image size

What remains unaddressed is what is powering this large mass-loss rate. One possibility is that the required energy is injected via the irradiation of the donor by the luminosity generated by the neutron star. For IGR J17062, the accretion luminosity is about 6 × 10³⁵ erg s⁻¹. Assuming isotropic emission, the fraction of radiation intercepted by the companion star is $f={({R}_{c}/2a)}^{2}=0.23 \%$, giving an irradiation rate of ${\dot{E}}_{\mathrm{irr}}=1.3\times {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. The induced mass-loss rate can be estimated as

$$\begin{eqnarray}&&\dot{M}=-\eta {\dot{E}}_{\mathrm{irr}}\displaystyle \frac{{R}_{c}}{{{GM}}_{c}},\end{eqnarray} \tag{ 11 }$$

so that η ≈ 3% of the impinging energy must be converted into the kinetic energy of the wind. In practice, this efficiency should be considered a lower bound, as the companion star is likely shadowed by the accretion disk, and any anisotropy in the accretion luminosity is unlikely to beam the emission into the orbital plane. Whether or not this efficiency is realistic is not clear, as the efficiency expected from theory will depend greatly on the energy spectrum of the wind-driving radiation and the geometry of the companion star. Typically assumed efficiencies range from 0.01%–10% (Tavani & London 1993; Justham et al. 2006; Ginzburg & Quataert 2020), so it seems that the efficiency required to drive a wind in IGR J17062 is at least plausible, if on the high side.

Finally, we note that Strohmayer et al. (2018) derive a system inclination of 19°–275 based on the assumption that the binary evolves on the timescale of gravitational wave emission. Clearly, this assumption cannot be correct, so that the constraint on the inclination does not hold either.

Consistent with previous spectral analyses of IGR J17062 (Degenaar et al. 2017b; Keek et al. 2017; van den Eijnden et al. 2018), we found that the NICER data were well described using a phenomenological model consisting of an absorbed blackbody and power law, with an additional Gaussian emission line centered at 1 keV. Resolving the energy spectra as a function of the pulsar rotational phase, we find a significant pulse modulation in the blackbody temperature and the power-law flux. Additionally, we found marginal evidence (2σ–3σ) for modulation in the blackbody normalization, the Gaussian line flux, and the power-law photon index.

One should be careful not to overinterpret the results from phase-resolved spectroscopy. The adopted spectral model is entirely phenomenological and likely does not fully reflect the physical processes underlying the observed emission. A more robust analysis of the phase-resolved spectroscopy requires detailed waveform modeling (see, e.g., Salmi et al. 2018), which is well beyond the scope of this work. With that caveat in mind, however, we note that the significant temperature oscillation of the blackbody component does support the interpretation that this thermal emission component is associated with the stellar surface (van den Eijnden et al. 2018; Hernández Santisteban et al. 2019).

A potentially quite interesting aspect of the phase-resolved spectroscopy is that there is some (weak) evidence that the Gaussian line flux is modulated at the stellar rotation rate. The precise origin of this emission line is not known, but it is likely with either Fe–L band emission, or possibly with a Ne ii line (Degenaar et al. 2013). High-resolution spectroscopy indicates the 1 keV feature can be resolved into several narrow lines (Degenaar et al. 2017b; van den Eijnden et al. 2018) and might be attributed to a collisionally ionized plasma, possibly from a shock in the accretion column. Alternatively, the line might be due to a blueshifted emission line in an outflow or from an accretion disk reflection feature (Degenaar et al. 2017b; van den Eijnden et al. 2018). A possible phase dependence of the line flux would offer interesting context to this debate. If the line originates in a shock, then one would indeed expect that the line is modulated by the stellar rotation, in which case the precise phase delay and waveform of the line flux should carry information on the precise location of that shock. Alternatively, if the line is associated with an outflow or a reflection feature, we can hypothesize that a phase dependence arises indirectly from the beamed pulsar emission sweeping periodically over the disk (or outflow), opening an opportunity for a tomographic reflection study of the accretion flow (see, e.g., Ingram et al. 2017). Given the long-term stability of its pulsations, the data quality of the energy-resolved pulsations can in principle be improved with further observations. Hence, we suggest that IGR J17062 is an especially interesting target for continued monitoring and detailed pulse waveform modeling.

Considering the pulse amplitudes measured in each of the separate observing epochs between 2016 and 2020, we found that the 0.4–6 keV pulse fraction of IGR J17062 ranges from 0.5% to 2.5% following a slow evolution with time. This evolution could be well described by a simple sinusoidal oscillation, yielding a period measurement of 1210 ± 40 days. Because we only observed a little over a single cycle, we could not determine if this apparent modulation is maintained over decades, nor if the oscillation is periodic or quasi-periodic.

The 2008 RXTE observation does not yield any useful constraints on the long-term pulse-amplitude evolution. While the 2–12 keV pulse amplitude of the RXTE observation has a much higher pulse fraction of 9.4% ± 1.1% (Strohmayer & Keek 2017), we found that this measurement is strongly dependent on photon energy (Figure 5). At 2 keV, the pulse fraction drops to only 6% ± 1%. If we compare this to the time-averaged pulse fraction found with NICER, we see that the measurements are consistent within the statistical uncertainty. What little tension still remains between NICER and RXTE observations could be due to the time dependence. If we extrapolate the periodicity observed over the past four years, then the RXTE observation coincides with a peak of the sinusoid and should be expected to yield a higher pulse fraction than the time-averaged NICER data. Alternatively, if the time modulation of the pulse fraction is not strictly periodic, then the residual tension could worsen, suggesting a long-term decay of the pulse fraction. Continued monitoring of the pulse fraction oscillation may be able to resolve this uncertainty by determining the stability of the period. In that context, we note that if the recent evolution holds, the next peak in pulse fraction should occur around 2021 September, with the subsequent minimum occurring around 2023 May.

The origin of the variable pulse amplitude is not clear at this time. Because the pulsations of an AMXP are powered by the variable accretion flow, it is generally expected that some of the accretion variability is visible in the pulse waveforms as well. Indeed, the sample of known AMXPs shows ample evidence for such a transfer mechanism. On rapid timescales (seconds), the stochastic and quasi-periodic noise of the accretion disk has been shown to modulate the pulse amplitude (Menna et al. 2003; Bult et al. 2017). Pulse amplitudes also routinely show variations on slower day-long timescales, which are generally understood to be stochastically driven, and are obviously correlated with variations in the source luminosity (Chou et al. 2008), spectroscopy (Kajava et al. 2011), or temporal variability (Bult & van der Klis 2015). The slowly oscillating pulse fraction of IGR J17062 appears to be of a very different nature. Its timescale is much longer than any pulse variability seen in other AMXPs, and notably, it does not correlate with changes in the pulse phase, source luminosity, or any other observable of the accretion system. In particular, the lack of correlation between the pulse fraction and the X-ray flux would suggest that the modulation is not driven by the mass accretion rate. Instead, it may originate in a geometric effect such as a periodic modulation of the apparent size and/or orientation of the stellar hot spot. For instance, it may be that the neutron star is slowly precessing (Zimmermann & Szedenits 1979; Alpar & Pines 1985; Jones & Andersson 2001; Link 2003), although the precession period in AMXPs has been suggested to be on the order of minutes to days (Chung et al. 2008) and presumably would impact the pulse phase as well.

Intriguingly, the only other process in IGR J17062 with a similar timescale is the duty cycle of its intermediate-duration X-ray bursts. Assuming no X-ray bursts were missed, the three bursts observed in 2012, 2015, and 2020 give an average burst recurrence time of 1460 days, only somewhat larger than the apparent cycle period observed in the pulse fraction. Both the 2015 and 2020 X-ray bursts occurred near the minima of the pulse fraction. For this relation to extend to the 2012 X-ray burst as well, however, the pulse-amplitude oscillation would have to be quasi-periodic, so continued monitoring should be able to establish if these phenomena are truly connected. If, however, they are indeed related, then we can speculate that the changes in the pulse fraction are tied to some physics associated with the deep neutron star crust, rather than the accretion flow.

A remaining open question is why, after two decades in outburst, the pulsations of IGR J17062 are still visible. The two other AMXPs that exhibited a multiyear outburst only showed pulsations for a limited time early in their outbursts, suggesting that perhaps the decaying pulsations are a feature of prolonged accretion episodes. We briefly discuss how these other two cases compare to IGR J17062.

The best-studied case is that of HETE J1900.1–2455 (HETE J1900). This source was first discovered in 2005 June (Vanderspek et al. 2005) and remained in outburst until late 2015 (Degenaar et al. 2017a). Throughout most of its ∼10 yr outburst, HETE J1900 maintained a high luminosity of about ∼4 × 10³⁶ erg s⁻¹ (Falanga et al. 2007; Papitto et al. 2013), with significant intensity variability throughout (Patruno & Wijnands 2017; Degenaar et al. 2017b). Its 377 Hz pulsations were persistently visible for about 22 days (Kaaret et al. 2006; Galloway et al. 2008), before becoming intermittent and subsequently disappearing entirely (Patruno 2012). This decay has been suggested to be due to magnetic field burial (Cumming 2008), which is supported by the apparent exponential decay in the accretion torque derived from the pulsar timing (Patruno 2012). A particularly strange aspect of HETE J1900 is that during the first two months of outburst, its pulse amplitudes appeared to be tied to the occurrences of X-ray bursts: following an X-ray burst, the pulse fraction would abruptly increase and then steadily decay over a ∼10 day timescale (Galloway et al. 2007).

In IGR J17062, the X-ray bursts might also be occurring when the pulse amplitude is at a minimum, but we did not observe a similar abrupt jump in the pulse fraction following the 2020 burst. The analysis of the NICER epoch containing the X-ray burst was complicated by the impact of the X-ray burst cooling tail, though. While the X-ray burst emission dominates, we would not expect to see accretion-powered pulsations (see, e.g., Watts 2012). The subsequent dip in source intensity suggests that accretion onto the neutron star is inhibited, in which case we might reasonably expect the pulsations to be likewise suppressed. Hence, we would expect that the pulse signal is significantly variable during this time. Unfortunately, however, we lack the sensitivity to test this theory. We find an upper limit on the pulse amplitude of 0.8%–1.1%, which is comparable to the pulse fraction we would expect from interpolating the long-term evolution. Hence, we cannot be certain if the pulsations are suppressed or if they are present at the expected rate. We can, however, rule out an abrupt factor of ≳2 jump in pulse fraction like the one observed in HETE J1900 (Galloway et al. 2007).

The second case is that of MAXI J0911–655 (MAXI J0911). First discovered in 2016 February (Serino et al. 2016), the outburst of this source is still ongoing. Observations collected during the first month in outburst found 340 Hz pulsations (Sanna et al. 2017b), but these pulsations were no longer visible in observations collected at later times (Bult et al. 2019). In addition to the longevity of its outburst, MAXI J0911 shares a couple of similarities with IGR J17062: it has a 44 minute ultracompact binary orbit (Sanna et al. 2017b), and it shows similarly rare but energetic intermediate-duration X-ray bursts (Nakajima et al. 2020; Bult et al. 2020b). Where the sources differ, though, is in their luminosity. Like HETE J1900, the luminosity of MAXI J0911 is estimated to be ∼5 × 10³⁶ erg s⁻¹ (Sanna et al. 2017b), an order of magnitude larger than the ∼6 × 10³⁵ erg s⁻¹ luminosity of IGR J17062. Hence, a simple explanation for the persistently visible pulsations of IGR J17062 might be rooted in the much lower luminosity. That is, the magnetic screening process may simply be ineffective at the much lower accretion rate of IGR J17062 (Cumming et al. 2001) or operate on much longer timescales.