Is the black-widow pulsar PSR J1555–2908 in a hierarchical triple system?

PSR J1555$-$2908, a neutron star spinning rapidly at $559$ Hz, resides in a tight $5.5$ hr binary system (Ray et al. 2022, submitted, hereafter Paper I) with a low-mass ($0.06\,M_{\odot}$) companion (Kennedy et al. 2022, submitted). This millisecond pulsar (MSP) was first detected in a Green Bank Telescope (GBT) pulsar survey which targeted steep-spectrum radio sources (Frail et al., 2018) within the localization region of Fermi-Large Area Telescope (LAT; Atwood et al., 2009) gamma-ray sources. After the radio discovery, gamma-ray pulsations were detected, allowing the timing measurement of the system parameters over the $12$-yr Fermi mission time span.

Timing analysis using the Large Area Telescope (LAT) data reveals variations of the spin frequency $f$ which are larger than is typical for MSPs (Kerr et al., 2015). Such variations are often seen in young gamma-ray pulsars, where they are labeled as “timing noise”. Timing noise is also present in MSP, but its weaker amplitude (Shannon & Cordes, 2010) is typically not detectable with the LAT (Kerr et al. 2022, submitted). Instead, the intrinsic rotational phase is generally well described by a quadratic function of time, i.e. a constant spin-down rate $\dot{f}$, with no detectable noise components above the measurement uncertainty. In contrast, for PSR J1555$-$2908 four additional terms in the rotational phase model are required.¹¹1In a gamma-ray timing analysis, timing noise is easily distinctive from the variations of the orbital period which are often seen in “spider” pulsars. While both variations happen on multi-year timescales, spin-frequency variations manifest as pulse phase drift, and variations of the short ($\lesssim 1$ d) orbital period smear the pulse profile out.

Such strong variations suggest an alternative explanation. For example, in the case of the millisecond pulsar PSR J1024$-$0719, long thought to be isolated, the measurement of additional higher-order spin-frequency derivatives led to the conclusion that the pulsar could be in a very-long-period ($>2$ kyr) orbit with a low-mass companion (Kaplan et al., 2016).

In this paper, we discuss a similar alternative to the timing noise explanation for PSR J1555$-$2908: that the variations arise from the presence of a third body in the system. This additional object is in a wide, multi-year orbit around the closely orbiting neutron star and its low-mass companion. An extreme example of such a system would be the hierarchical-triple-system pulsar PSR J0337$+$1715 (Ransom et al., 2014).

With the currently available data, this model describes the pulsar as well as the timing-noise model, but provides a simple and clear physical explanation for the spin frequency variations.

To precisely track the rotational phase, the photon arrival times need to be corrected for the line-of-sight motion of the pulsar, if it is in an orbit with one or more companions. For a simple hierarchical triple system (HTS), we assume that the gravitational interaction between the two companions can be neglected, and the third body orbits the center of mass of the inner binary system. The photon times, corrected for the pulsar’s motion around the barycenter of the triple system, $t_{\rm tb}$, can be expressed as a function of the photon’s emission time, $t_{\rm em}$, as

Here, $x_{\rm i}$ and $x_{\rm o}$ are the times light needs to travel the semimajor axis of the pulsar’s orbit due to the inner (i) and outer (o) companion projected onto the line of sight. $\Delta_{\rm i}$ and $\Delta_{\rm o}$ are dimensionless functions that describe the respective modulations depending on the orbital phases, $\Omega_{\rm i}t_{\rm em}$ and $\Omega_{\rm o}t_{\rm em}$. The orbital angular frequencies $\Omega_{\rm i}$ and $\Omega_{\rm o}$ are defined by the orbital periods $P_{\rm i}$ and $P_{\rm o}$ via $\Omega_{\rm i}=2\pi/P_{\rm i}$ and $\Omega_{\rm o}=2\pi/P_{\rm o}$.

Following Edwards et al. (2006), we can make $x_{\rm i}\Delta_{\rm i}(\Omega_{\rm i}t_{\rm em})+x_{\rm o}\Delta_{\rm o}(\Omega_{\rm o}t_{\rm em})$ a function of $t_{\rm tb}$ by Taylor expansion. To second order (i.e. including terms of order $(x_{\rm i}\Omega_{\rm i})^{2}$, $(x_{\rm o}\Omega_{\rm o})^{2}$, $x_{\rm i}x_{\rm o}\Omega_{\rm i}^{2}$, and $x_{\rm i}x_{\rm o}\Omega_{\rm o}^{2}$), we find

where $\Delta^{\prime}$ and $\Delta^{\prime\prime}$ denote the first and second derivative of $\Delta$ with respect to $\Omega t$. For wide, multi-year outer orbits of low-mass companions, $x_{\rm o}\Omega_{\rm o}\ll x_{\rm i}\Omega_{\rm i}$, so we neglect terms multiplied by this quantity, i.e. $\Delta_{\rm o}^{\prime}=0$ and $\Delta_{\rm o}^{\prime\prime}=0$.

In the case of a circular orbit, $\Delta$ takes the simple form of a sinusoid, $\Delta(t)=\sin(\Omega(t-T_{\rm asc}))$, with time of ascending node $T_{\rm asc}$. Hence, $\Delta^{\prime}(t)=\cos(\Omega(t-T_{\rm asc}))$ and $\Delta^{\prime\prime}(t)=-\sin(\Omega(t-T_{\rm asc}))$. Eq. (2) can also be applied to the widely used small-eccentricity-orbit “ELL1 model” (Lange et al., 2001), or the larger-eccentricity models presented in Nieder et al. (2020).

For the outer companion’s orbit of PSR J1555$-$2908, we use the full “BT model”, $\Delta(t)=\sin\omega(\cos E-e)+\sqrt{1-e^{2}}\cos\omega\sin E$ with eccentricity $e$, and angle of periastron $\omega$ (Blandford & Teukolsky, 1976). $E$ is defined via $E-e\sin E=\Omega_{\rm b}(t-T_{0})$ and is computed numerically. Instead of using the parameter set $\{T_{0},e,\omega\}$, which is degenerate for small eccentricities, we use $T_{\rm asc}=T_{0}-\omega/\Omega_{\rm b}$, $\epsilon_{1}=e\sin\omega$, and $\epsilon_{2}=e\cos\omega$ (Lange et al., 2001).

In Section 3, we show that PSR J1555$-$2908 can be well described as an HTS with the inner orbit being circular and the outer orbit being eccentric.

For a timing analysis, it is helpful to have a good starting point in the multi-dimensional parameter space. The parameters describing the potential outer orbit can be roughly estimated from the timing solution presented in Paper I, which requires four additional frequency derivatives (“FX model”).

First, we need to derive how the observed frequency changes in the case of an outer circular orbit. Assuming only a non-zero spin frequency $f$ and a non-zero first time-derivative $\dot{f}$, the additional evolution of the spin frequency over time can be written in two ways:

Here, we denote the $n$th frequency derivative as $f^{(n)}$, and use an asterisk to indicate parameters measured with the FX model, and $\Delta f_{*}=f-f_{*}$ as well as $\Delta\dot{f}_{*}=\dot{f}-\dot{f}_{*}$. $t_{0}$ denotes the chosen reference time.

Second, since all parameters in Eq. (3) are measured, we fit Eq. (4) to Eq. (3) over the valid time span of the pulsar timing solution. To do this, we used curve_fit from scipy. This gives initial values for the five free parameters $\{\Delta f_{*},\Delta\dot{f}_{*},x_{\rm o},\Omega_{\rm o},T_{\rm asc,o}\}$ from Eq. (4).

From the fitted values, we find the inner binary barycenter’s potential movement would have a radius of a few tens of kilometers with an orbital period of $\sim 10$ yr, indicating a low-mass companion in a wide orbit.

In our timing analysis, we use the LAT data prepared for the timing in Paper I. These data included P8R3 SOURCE-class photons (Atwood et al., 2013; Bruel et al., 2018) detected by LAT between 2008 August 03 and 2020 August 05 within a $5\deg$ region of interest (RoI) around the pulsar position, with energies greater than $100\,$MeV, and with a maximum zenith angle of $90\deg$.

The timing analysis is done twice, very similarly to the approach in Paper I. Firstly, we refit all the parameters of the FX model, i.e. we did not fix any parameters to the radio timing solution. Secondly, we amended the timing code with the HTS model given in Equations (1) and (2). To judge the timing results, we used the $H$ statistic (de Jager et al., 1989; Kerr, 2011) with the photon probability weights optimized following Bruel (2019), the likelihood $\mathcal{L}$ (Abdo et al., 2013), and the Bayesian Information Criterion (BIC; Schwarz, 1978). The latter is based on the likelihood but with a penalty factor proportional to the degrees of freedom to penalize higher-dimensional models.

In the timing analysis, we tried to find the parameters that maximize $\mathcal{L}$ for a given model. In the case of the HTS model, we tested the circular and eccentric options, where the latter requires two additional parameters. The eccentric orbit model is favored by all three statistics. For the FX model, we added spin-frequency derivatives until the test statistics $H$ and $\mathcal{L}$ stopped improving significantly, which is in turn disfavored by the Bayesian Information Criterion (BIC). Here, adding six additional derivatives in principle leads to $H$ and $\mathcal{L}$ values of a similar size to the HTS model. However, the last model favored by the BIC is the one with four additional frequency derivatives. The final timing solutions for the FX and HTS models are reported in Table 1.

In the timing analysis with the HTS model, the highest test statistics are found for an eccentric ($e_{\rm o}\sim 0.27$) outer orbit of $\sim 4{,}500$ d. With the LAT data span being only slightly shorter, the spin parameters $f$ and $\dot{f}$ are still correlated with the outer-orbit Keplerian parameters $P_{\rm o}$, $x_{\rm o}$ and $T_{\rm asc,o}$. However, this covariance is small enough to not pose a problem for the timing analysis. Note, that the timing analysis does not give any information on the inclination between the orbits as only the line-of-sight motion can be measured.

The highest test statistics, $H$ and $\mathcal{L}$, seem to favor the HTS model over the FX model (see Table 1). Since each model has a different number of free parameters, we invoke the BIC. The BIC’s penalty for the additional parameter removes most of the advantage the HTS model had. However, the BIC still slightly favors the HTS model over the FX model ($\Delta\textrm{BIC}=-1.6$). Our final timing solutions over $12$ years of LAT data using the FX model and the HTS model are shown in Table 1.

In order to more directly compare the timing noise to that of other MSPs, we have analyzed it by following the methods of Lentati et al. (2014), adapted to gamma-ray data as in Kerr et al. (2022, submitted). In this approach, the timing noise is assumed to originate in a stationary process, thus having a unique description with a power spectral density which we assume to follow a power law, $P(f)=P_{\mathrm{tn}}(f\times\mathrm{yr})^{-\Gamma_{\mathrm{tn}}}$. The timing model parameters, timing noise parameters $\Gamma$ and $P_{\mathrm{tn}}$, and the Fourier coefficients of the timing noise process are all determined jointly by maximizing the likelihood, which includes a Gaussian component that constrains the Fourier coefficients to follow the specified power law. We estimate uncertainties on the timing noise parameters from the curvature of the likelihood at its maximum value, and obtain $\log_{10}(P_{\mathrm{tn}}\times\mathrm{s^{-2}\,yr})=-13.3\pm 1.3$ and spectral index $\Gamma_{\mathrm{tn}}=6.0\pm 1.4$.

The spin-frequency variations of PSR J1555$-$2908 (see also Paper I) are unusually large in amplitude for the timing noise of an MSP (Kerr et al., 2015). Indeed, Figure 1 shows the timing noise spectral properties for two samples of MSPs as measured by the pulsar timing arrays NANOGrav (Arzoumanian et al., 2020) and the PPTA (Goncharov et al., 2021), along with our measured values for PSR J1555$-$2908. It is possible this collection of timing noise measurements is biased by selection for inclusion in a Pulsar Timing Array (PTA). However, the driving concern for inclusion is radio brightness, and often timing noise emerges only after many years of observation. Consequently, we expect such a bias to be modest. The abscissa in Figure 1 is defined in such a way to indicate the timing noise at frequency $f=0.1\,\mathrm{yr}^{-1}$, i.e. a time scale comparable to the length of the data set. PSR J1555$-$2908 is a clear outlier, with more than $10$ times the timing noise amplitude of the next strongest pulsar, J1824$-$2452A in the M28 globular cluster, whose timing noise presumably includes contributions from the cluster potential and close passages by other cluster members. Another noteworthy example is PSR J1939$+$2134, an isolated MSP considered to be a “noisy” timer for which Shannon et al. (2013) invoked an asteroid belt surrounding the pulsar to explain the observed timing variations, but which has a timing noise amplitude that is $1000$ times weaker than that of PSR J1555$-$2908.

In this work we have presented an alternative explanation, in which another companion orbits the inner binary system, a hierarchical triple system. To account for the additional body in the system, we have developed a simple, but sufficient rotational phase model. In a timing analysis with current data, the results with both models are indistinguishable. Within the timing range the maximum phase offset between both models is $<1\%$ of one rotation.

It is predictable that the two models would not differ much. The variations are measurable with the LAT data, but the peak-to-peak amplitude of the phase offset is still only $17\%$ of one rotation. Furthermore, the potential outer companion would have completed roughly one orbit during the Fermi mission which can be well approximated with the six spin parameters of the FX model.

In the following, we will make estimates of the companion mass and the orbital radius, but want to note that these directly depend on the outer orbital period $P_{\rm o}$. Due to the strong correlations between the orbital parameters (see Section 3), these estimates should be treated with caution and will only give an order of magnitude.

The binary mass function may be used to estimate the mass of the outer companion (see Fig. 2). For common pulsar masses $1.4\,M_{\odot}<M_{\rm p}<2.0\,M_{\odot}$, an inner companion mass $M_{\rm i}=0.06\,M_{\odot}$ (Kennedy et al. 2022, submitted), and inclination angles $~{}30\arcdeg-90\arcdeg$, the outer companion would have roughly a Mercury-like mass $0.02\,M_{\oplus}\lesssim M_{\rm o}\lesssim 0.06\,M_{\oplus}$. To date, this companion would still be listed among the four lowest-mass extrasolar planets²²2http://exoplanet.eu/catalog/.

A comparably low-mass ($\sim 0.015\,M_{\oplus}$) companion has also been found for a planet orbiting PSR B1257$+$12 (Wolszczan, 1994), the pulsar thanks to which, $30$ years ago, the first extrasolar planetary companions have been identified (Wolszczan & Frail, 1992). However, unlike the periods of the planets orbiting PSR B1257$+$12 with tens to hundreds of days, the probable companion of PSR J1555$-$2908 takes considerably longer with $P_{\rm o}\sim 4{,}500$ days.

The putative triple system of PSR J1555$-$2908 is extremely hierarchical. The pulsar’s mass is $\sim 28$ times larger than the inner companion’s mass, and even $\sim 10^{7}$ times larger than the outer companion’s mass. The ratio between outer and inner orbital period is $P_{\rm o}/P_{\rm i}\approx 20{,}000$. Using Kepler’s third equation, the masses from optical modeling, and the timing parameters one can estimate the semimajor axis of the inner companion’s orbit as $\sim 4.2$ lt-s and the outer companion’s orbit as $\sim 3{,}200$ lt-s. Systems with such extreme hierarchical properties are typically assumed to be stable on long timescales (Georgakarakos, 2008).

Other HTSs which include a pulsar are known. PSR J0337$+$1715 is in an orbit with two white dwarf companions (Ransom et al., 2014). Carefully studying this system led to one of the most stringent tests of General Relativity’s predicted universality of free fall (Archibald et al., 2018; Voisin et al., 2020). A binary system consisting of the MSP PSR B1620$-$20 and a white dwarf is orbited by a planetary companion of $\sim 2.5\,M_{\rm Jupiter}$ (Thorsett et al., 1993; Sigurdsson et al., 2003).

An open question is how such a system could have formed. A pulsar planet can form in various ways (see, e.g., Podsiadlowski, 1993). Since PSR J1555$-$2908 is not in a dense Globular Cluster, it seems unlikely that the triple was formed in an exchange interaction between systems as is suggested for PSR B1620$-$20 (Sigurdsson et al., 2003). Another option could be that the inner binary captured a free-floating planet (see, e.g., Miret-Roig et al., 2021). On the other hand, the planet might have formed from a circumbinary debris disk (see, e.g., Pelupessy & Portegies Zwart, 2013), in which case it would be expected that the orbits are coplanar.

Only to emphasize the power of pulsar timing with LAT data, we want to note that if PSR J1555$-$2908 has an additional planetary-mass companion, then it was discovered because over the course of $100$ billion pulsar rotations, the phase drifted by $17\%$ of one single rotation.

So, is PSR J1555$-$2908 in a hierarchical triple system? We do not know, yet. With current data, both of the timing models that are presented here track the pulsar’s rotation equally well. However, the ongoing Fermi mission will continue to collect data, and eventually the nature of this system will be resolved.