© The Authors 2025

1. Introduction

Being extremely accurate clocks (Manchester 2017), millisecond pulsars (MSPs) undergo particular scrutiny for any timing irregularity that may be the signature of a wide variety of phenomena, including studies of stellar evolution (e.g. Tauris & van den Heuvel 2023; Tauris et al. 2017), the discovery of the first extra-solar planets (Wolszczan & Frail 1992), the ongoing pursuit of a low-frequency gravitational-wave background (e.g. Reardon et al. 2021; Alam et al. 2021; Chen et al. 2021), constraining the neutron-star equation of state (e.g. Fonseca et al. 2021; Riley et al. 2021), or tests of gravity theories (e.g. Freire & Wex 2024; Kramer et al. 2021). One of the most prominent systems for the study of stellar evolution and especially for tests of gravity theories is PSR J0337+1715 (J0337 hereafter). This pulsar was discovered in the GBT drift-scan survey (Boyles et al. 2013; Lynch et al. 2013), and the triple nature of the system was confirmed by Ransom et al. (2014). The pulsar follows a ∼1.6-day orbit with the inner ∼0.2 M_⊙ helium white dwarf, which is detectable at optical wavelengths (Kaplan et al. 2014). Together they form what we call the inner binary. The centre of mass of the inner binary can be seen as orbiting with the outer ∼0.4 M_⊙ white dwarf in ∼327 days. Interestingly, the orbits have low eccentricities and are nearly co-planar, which has proved fundamental for understanding the evolution of the system (Tauris & van den Heuvel 2014).

Altogether, the triple system extends over ∼1 AU, which makes it uniquely compact compared to other multiple systems. The first pulsar that was found in a triple system was PSR B1620−26 (Backer et al. 1993; Thorsett et al. 1999). The system is composed of a likely gas giant companion orbiting a pulsar white dwarf binary in a hierarchical orbit (Sigurdsson et al. 2003; Sigurdsson & Thorsett 2005). To our knowledge, there are three other systems suspected to be hierarchical triple systems: ZTF J1406+1222, which is a black-widow pulsar candidate that has recently been suggested to be associated with a distant low-mass stellar companion in a ∼12 000-year orbit (Burdge et al. 2022) extending over ∼600 AU; the redback pulsar PSR 2043+1711, which has a suspected distant stellar companion (Donlon et al. 2024); and the pulsar binary PSR J1618−3921, which shows anomalous orbital-period variations and a spin-period second derivative compatible with the presence of a third object in the system(Grunthal et al. 2024). There is yet another system, PSR J1903+0327 (Champion et al. 2008), that is thought to have formed in a triple system that later became unstable (Freire et al. 2011; Portegies Zwart et al. 2011; Pijloo et al. 2012) and is currently a unique type of binary system.

The presence of an MSP in such a relatively compact triple stellar system allows for a unique test of the universality of free fall (UFF) involving a strongly self-gravitating object (Freire et al. 2012; Ransom et al. 2014; Shao 2016; Archibald et al. 2018; Voisin et al. 2020b, 2022). We now discuss the reasons for this and why it is important (for a more detailed explanation see Voisin et al. 2020b).

The UFF, extended to self-gravitating objects, is the gravitational weak equivalence principle (GWEP; Will 2014). Together with the local Lorentz and position invariances of gravity, it constitutes the strong equivalence principle (SEP). Testing the SEP is of great importance because of all valid theories of gravity, only general relativity (GR) fully embodies it¹; all others predict some type of SEP violation. This means that a violation of the SEP is perhaps the most promising avenue for finding new physics beyond GR. This also means that searching for SEP violation will either rule out GR (in case of a detection) or alternative theories of gravity (in case of a non-detection).

In alternative theories of gravity, such as scalar-tensor theories (Damour & Esposito-Farèse 1993; Damour & Esposito-Farese 1996), violation of GWEP takes the form of an effective body-dependent gravitational constant G_ab = G_ba = G(1 + Δ_ab), where G is the Newtonian gravitational constant and Δ_ab is a violation parameter, the definition of which is theory dependent. In weakly self-gravitating bodies this reduces to the usual discrepancy between inertial and gravitational masses. Assuming WEP applies, the ratio of the two masses is then proportional to the fractional gravitational binding energy, where the proportionality coefficient is the so-called Nordtvedt parameter (Nordtvedt 1968). One can see that in a binary system orbiting according to Newton’s equations of motion, substituting G for G_ab only amounts to re-scaling masses, and with those being unknown, it does not allow for any constraint on the SEP (which is no longer true with post-Newtonian equations of motion). In a triple system, however, there is no such degeneracy. That is why the WEP has been tested using two masses of various compositions falling in the gravitational field of a third one, as in the case of the MICROSCOPE experiment, which had two masses in Earth orbit for one year (Touboul et al. 2022).

The GWEP was first tested with the Lunar Laser Ranging (LLR) experiment, where the triple system of self-gravitating objects is the Earth-Moon-Sun (⊕ −°− ⊙) system. The extremely accurate determination of the Moon’s distance provided by LLR resulted in a tight limit on the Nordtvedt effect (Nordtvedt 1968) and a resulting limit of Δ_⊕⊙ − Δ_°⊙ = (3.0 ± 5.0)×10⁻¹⁴ (Hofmann & Müller 2018). However accurate, this result provides a weak limit to the SEP (in particular the Nordtvedt parameter) because of the small binding energies of the Earth and Moon. Furthermore, such weak-field limits do not constrain the strong-field regime because of additional theory-dependent effects that may arise, such as spontaneous scalarisation (Damour & Esposito-Farèse 1993; Damour & Esposito-Farese 1996).

For this reason, experiments with pulsars as the compact object are important. Prior to the discovery of J0337, the best that could be done was the Damour-Schäffer test (Damour & Schäfer 1991), where the triple system is effectively a near-circular pulsar – white dwarf binary evolving in the gravitational potential of the Galaxy. A violation of GWEP would polarise the orbital eccentricity of the binary towards the Galactic centre. Using one of the most accurately timed pulsars, PSR J1713+0747, Zhu et al. (2019) obtained |Δ_pG − Δ_cG|< 2 × 10⁻³ at a 95% confidence level (95% C.L.) from the lack of variation of the orbital eccentricity of the system. This type of test has one main limitation: the weak gravitational acceleration of the Galaxy. After the discovery of J0337, it became clear that the system is an extremely sensitive probe of GWEP violations given the (locally) much stronger field of the outer white dwarf star compared to the field of the Galaxy.

As was first shown in Ransom et al. (2014) and Archibald et al. (2018), the aforementioned compactness of the J0337 system makes it necessary to numerically integrate the three-body equations of motion at first post-Newtonian order to model the timing data. This represents a large computational cost compared to the usual pulsar-timing modelling, which relies on analytical expressions for binary systems. To perform these numerical integrations, Voisin et al. (2020b) independently developed a numerical timing model (Nutimo) (Voisin et al. 2020b) that supports the strong-field generalisations of the parametrised post-Newtonian equations of motion (for instance Damour & Taylor 1992 and Appendix A of Voisin et al. 2020b), which is necessary for a test of alternative theories of gravity, particularly the UFF. J0337’s orbital motion is weakly relativistic, and as such, it is not very sensitive to post-Newtonian effects (Damour & Taylor 1992), which are the key to test GR in relativistic pulsar binaries (Kramer et al. 2021); this is not a problem because any violation of the GWEP would appear at the Newtonian level in the equations of motion.

Using independent datasets and independent analysis, Archibald et al. (2018) and Voisin et al. (2020b) showed that this system constrains deviations from GWEP to |Δ|< 2.6 ⋅ 10⁻⁶ and |Δ|< 2.1 ⋅ 10⁻⁶, respectively, at a 95% C.L., where Δ denotes a violation of GWEP between the pulsar and any of the companions. In both cases, these results represent an improvement of nearly three orders of magnitude compared to the best previously available test of GWEP with a neutron star, using the above-mentioned PSR 1713+0747 (Zhu et al. 2015). This demonstrates how relevant J0337 is for this kind of work. A notable difference between the two approaches is the fact that the limit of Archibald et al. (2018) is mostly based on limits on systematics, while the limit of Voisin et al. (2020b) is mostly statistical. The latter was additionally able to constrain the longitudes of ascending nodes of the system.

Motivation of this work. Here, we continue the work of Voisin et al. (2020b) with an extended timing dataset of J0337 from the Nançay radiotelescope, now spanning approximately eight years (against five previously). The main motivation for the continued timing of this system is to improve the test of the GWEP provided by J0337.

A second motivation for this paper, which is linked with the first one, is to investigate an unmodelled low-frequency signal in the timing of J0337. This signal was already present in Voisin et al. (2020b), but it was interpreted as a possible red-noise process intrinsic to the pulsar emission mechanism (e.g. Shannon & Cordes 2010; Lyne et al. 2010) and left unmodelled. Instead, in Voisin et al. (2020b) uncertainties were widened in order to accommodate this effect and produce a conservative estimate of the posterior distribution functions, particularly that of the Δ parameter. However, in the more recent data discussed in this work, its amplitude has been increased to ∼4 μs, which is to be compared to the ∼2 μs uncertainty on individual times of arrival. This amplitude is such that it systematically limits the GWEP test with J0337. Thus, we aim to model it and evaluate the nature of its cause.

To do this, we considered two types of models for this low-frequency component: (1) red timing noise (chromatic or achromatic) and (2) the presence of a planet in a wide hierarchical orbit around the triple stellar system. The latter option is motivated by the unusually large amplitude of the residual signal compared to other known red-noise processes in millisecond pulsars (see Sect. 5). Precedents of companions initially associated with strong timing noise include the already mentioned planet around the pulsar binary PSR B1620−26 (Backer et al. 1993; Thorsett et al. 1999) as well as the likely distant stellar companion to the isolated pulsar PSR J1024−0719 (Bassa et al. 2016; Kaplan et al. 2016). To these can be added the proposal that the noise of PSR J1939+2134 (PSR B1937+21) may be caused by an asteroid belt Shannon et al. (2013). In any case, the apparent periodicity of the signal is comparable with the eight-year observation span, which will render our conclusions concerning its cause necessarily preliminary.

A third motivation for this paper is obtaining a more detailed study of the evolution of the system. Notably, this is also important for evaluating whether the presence of a planet is possible.

Since the discovery of PSR B1257+12 (Wolszczan & Frail 1992) and its system of two terrestrial-mass and one moon-mass planets (Konacki & Wolszczan 2003), only five more pulsars have been found to host planets (Nitu et al. 2022; Vleeschower et al. 2024, and references therein), four of which are dense Jupiter-mass ‘diamond planets’ that are likely remnant white dwarf cores (Bailes et al. 2011). The fifth is the planet in the PSR B1620−26 system (Backer et al. 1993; Thorsett et al. 1999; Sigurdsson et al. 2003). In addition are the already mentioned candidate companion to PSR J1555−2908 Nieder et al. (2022) and the possible asteroid belt around PSR J1939+2134 (PSR B1937+21) Shannon et al. (2013). Surveys have put limits on the planet population, particularly Kerr et al. (2015) around young pulsars and Behrens et al. (2020) for millisecond pulsars, both of which detected no planets.

More recently Nitu et al. (2022) studied a broader sample of 800 pulsars and found 15 candidates, most of which are considered likely to be quasi-periodic noise of magnetospheric origin. In some of these cases simultaneous pulse-profile variations can be observed, clearly pointing at that explanation.

These examples show the particular difficulty involved in detecting small planets. However, in the case of J0337, a clearer signature can be expected due to the specific mutual interactions within a four-body system.

In the following sections of this paper, we present the extended dataset in Sect. 2, describe the models in Sect. 3, and present the posterior inference and fitting results in Sect. 4. We end with a discussion of the consequences of the two main hypotheses, red noise or planet, in Sect. 5.

2. Observations

The pulsar J0337+1715 has been regularly observed since July 2013 every two or three days with the Nançay radio telescope (NRT) using its L-band receiver at a central frequency of 1484 MHz. The NRT is a meridian Kraus design collector equivalent to a 94-meter dish able to conduct ∼1 hour observations on any given source within its declination range each day. The dual linear polarisation signals are sent to the Nançay Ultimate Pulsar Processing Instrument (NUPPI, Desvignes et al. 2011; Cognard et al. 2013), an instrument that is able to coherently dedisperse a total bandwidth of 512 MHz. In this work, we use a dataset of 12474 20-mins integrated times of arrival (ToA) divided in four 128 MHz bands observed between MJD 56492 and MJD 59480 (July 2013 and September 2021)². As in the previous analysis presented in Voisin et al. (2020b) where details can be found, the times of arrival are estimated using the pat tool from the PSRCHIVE software library (Hotan et al. 2004), with the Fourier domain with Markov chain Monte Carlo (FDM) method. We note that two more years of data were added for this analysis.

At MJD 58631, new cooled pre-amplifiers were installed at the telescope, boosting the sensibility in the upper part of the band (∼1550–1730 MHz) and slightly changing the overall pulse shape (the exact pulse shape is smoothly frequency dependent). Around MJD 58780, we started a better polarisation calibration observing a slow, strong and polarised pulsar over one hour and the receiver rotating for ∼180 degrees (Guillemot et al. 2023). Thus the quality of the data determination in the range MJD 58631–58780 has been altered due to inappropriate calibration. In practice, this leads to significantly larger ToA uncertainties during that ∼150 days period, reflecting that small discrepancy between the template and the distorted daily profiles. As a result, although we do not expect the timing analysis to be significantly biased (thanks to wider uncertainties), we preferred to conservatively excise the range MJD 58631–58780 from our dataset.

3. Models

The present analysis is based on the numerical timing model of Voisin et al. (2020b). In the remainder of this section, we focus on two types of extension of the model aimed at explaining the residual low-frequency signal (see footnote²). The first extension assumes that the signal is due to a generic achromatic red-noise process modelled by a Fourier decomposition with a power-law spectrum, and an optional chromatic component caused by non-linear variations of dispersion-measure (DM) over time. The second extension assumes that the signal is caused by a small planet in a hierarchical orbit with the triple system. Parameters associated with the triple system are defined as in Voisin et al. (2020b), but nonetheless reviewed in Sect. 3.2 such that this paper is self-contained.

All models include astrometric and pulsar spin parameters that were described in Voisin (2017); Voisin et al. (2020b). Astrometric parameters are right ascension, declination and distance as well as the associated proper motion. Intrinsic parameters are the spin frequency and its derivative re-scaled to absorb the linear and quadratic effects of the Einstein and Shklovskii delays (Voisin et al. 2020b). Except in the dedicated chromatic model (see below), dispersion measure is fitted using two parameters, a constant value and a linear drift.

3.1. Red noise

Red noise in pulsar timing can be due to several causes. If it is chromatic, then the most common cause are variations of dispersion measure, which are a propagation effect caused by fluctuations of the column density of electrons in the interstellar medium. Such variations can be disentangled from achromatic signals if the observation is spectrally resolved, which is the case in the present work. Other sources of chromatic noise are scattering and band noise, which we do not consider in this work. Achromatic red noise, on the other hand, is intrinsic to the system (e.g. Shannon & Cordes 2010), and has been proposed to be caused by magnetospheric variations Lyne et al. (2010); Tsang & Gourgouliatos (2013), changes in the neutron star interior (e.g. Melatos & Link 2014, and references therein), or by an asteroid belt around the pulsar (Shannon et al. 2013).

The latter case illustrates the difficulty there can be in disentangling possible orbital motion caused by low-mass objects from other sources of noise. A recent example is the proposed presence of a small planet orbiting the black-widow system hosting PSR J1555−2908 (Nieder et al. 2022) in order to explain an unusually strong red noise signal.

3.1.1. Achromatic red noise

In recent years, red noise has been modelled for an increasing number of millisecond pulsars in the context of pulsar timing arrays (e.g. Chalumeau et al. 2022; Alam et al. 2021; Goncharov et al. 2021). Achromatic red noise is usually modelled as a truncated Fourier series the coefficients of which follow a Gaussian stochastic process with a power-law power-spectrum density. In order to select the best noise model, the posterior distribution function is marginalised over the deterministic part (the rest of the timing model) common to all models and selection is based on the comparison of their evidence using a Bayes factor (see, e.g. Chalumeau et al. 2022).

We faced several difficulties in the present work: (i) We wanted to compare the red-noise model to other models that, involving a planet, are deterministic. (ii) We could not perform analytical marginalisations over all the common parameters due to the numerical nature of the model. (iii) Evidence computation is computationally expensive. We note that point (iii) is largely a consequence of point (ii). On the other hand, given the relatively limited observation span, only a few Fourier components are actually necessary, which makes it possible to describe red noise with a deterministic model.

Thus, we added to the timing formula a truncated Fourier series the amplitudes of which follow a power law. This is equivalent to the average spectrum produced by a stochastic Gaussian process as described above. Formally, it can be written as

$$\begin{array}{c}\hfill F(t_a)=\sum _{k=1}^n\frac{A}{k^{\gamma }}sin(2k\pi \nu t_a+{\varphi }_k),\end{array}$$(1)

where A is the amplitude of the Fourier component at the fundamental frequency ν, γ is the power-law index, ϕ_k are phases at each harmonic k, and t_a are the times of arrival in the solar-system-barycentre frame. We note that the difference introduced by using times of arrival instead of times of emission is negligible given the time scale and amplitude of the signal.

We considered models from two to five Fourier components that we call ‘PLn’ where n ∈ {2, 3, 4, 5}. There are 3 + n fitted parameters, which are ν, A, γ, ϕ_{k ∈ [1, n]}. We remark that, contrary to the stochastic process approach, we fit for the fundamental frequency ν. Indeed, under the prior assumption that the Fourier spectrum is given by a power-law, the stochastic approach evaluates the probability of the signal being described by a complete Fourier series and therefore the fundamental frequency can be fixed to the inverse of the data span. Equation (1) is only a truncated series that therefore has additional freedom. Additionally, assuming the signal results from a stochastic process, here we fit a particular realisation of it that is unlikely to follow an exact power-law spectrum, especially one that is truncated. Fitting for the fundamental frequency may partially compensate for this.

3.1.2. Chromatic contribution: DMX

In order to test if the measured signal is due to variations in the interstellar medium, we tested a model combining an achromatic red-noise signal as described above with a variable dispersion measure (DM). The observation span was split in ten evenly distributed intervals each with a different DM_X value of DM, where X runs from one to ten. In this way, each interval covers ≃299 days, which is much smaller than the residual signal timescale of ∼3000 days and can capture it reasonably well if the signal turns out to be chromatic. On the other hand, this time interval is comparable to the orbital period of the outer binary (P_O ≃ 327 days), which avoids any correlations by averaging any effect due to the outer-binary motion over the interval.

3.2. A planet in a hierarchical orbit with the triple stellar system

A smooth, slow and quasi-sinusoidal signal is expected from a planet of mass m_π orbiting the triple system in a hierarchical orbit, that is with a period P_Π ≫ P_O with P_O the orbital period of the outer white dwarf, and a mass sufficiently small to ensure that its gravitational field is negligible compared to the fields within the triple system. In the present case, it is sufficient to consider planets of mass m_π ≪ M_⊙. The only measurable effect given the amplitude of the signal (a few μs) is the so-called Rømer delay induced by the planet, that is the variation of the distance between the pulsar and the observer induced by the presence of the planet.

One of the key difficulties is to disentangle the planet signal from red noise. Indeed, for relatively wide orbits (a few years) the Rømer delay of a low-mass companion is purely sinusoidal (with a first harmonic in case of an eccentric orbit), which can easily be absorbed into red-noise provided the signal is weak enough (due to low mass) and the number of observed orbital revolution is few, which for long wide orbits is usually the case. In the case of J0337, however, the more complex orbital configuration can in principle lift this degeneracy, as we show below.

In practice, we have implemented in Nutimo the possibility to add extra bodies the orbits of which are integrated numerically along with the orbits of the pulsar and the two white dwarfs. As in Voisin et al. (2020b), we have validated the accuracy of the integration in order to obtain a numerical timing accuracy at the level of a few nanoseconds.

3.2.1. Parametrisation

We extend here the orbital hierarchical parametrisation of the triple system (Voisin et al. 2020b; Archibald et al. 2018) in order to include a fourth body. In doing so, we explain how this parametrisation is related to Jacobi coordinates.

Let (R,V)_{k∈{p,i,o,π}} be the positions and velocities of the pulsar (p), inner white dwarf (i), outer white dwarf (o), and planet (π) relative to the inertial reference frame associated with their centre of mass. We call ‘inner binary’ the set I = {p,i} with centre of mass b, ‘outer binary’ the set O = {b,o} with centre of mass (t), which is also the centre of mass of the triple system, and we call ‘planetary binary’ the set Π = {t,π}, the centre of mass of which is also centre of mass of the whole system.

One can decompose the position of the pulsar as

$$\begin{array}{c}\hfill {\mathit{R}}_{\mathrm{p}}={\mathit{R}}_{\mathrm{p}/\mathrm{b}}+{\mathit{R}}_{\mathrm{b}/\mathrm{t}}+{\mathit{R}}_{\mathrm{t}},\end{array}$$(2)

and similarly its velocity V_p. In Voisin et al. (2020b), only the first two terms were present on the right-hand side. They describe the motion of the pulsar relative to the inner-binary centre of mass, and of the inner binary with respect to the centre of mass of the system. Here, we add a third term that describes the motion of the triple system with respect to the centre of mass of the whole system, which includes a planet. Thanks to the hierarchy of the system each term follows an approximately Keplerian motion (see also Appendix B), which makes a parametrisation of positions and velocities in terms of osculating orbital elements relevant. In addition, mutual interactions as well as Shapiro and Einstein delays allow one to constrain the four masses. This leads to the following correspondence between the state vectors and masses on one side, and orbital elements on the other side,

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{(\mathit{R},\mathit{V})}_{\mathrm{p}/\mathrm{b}},\\ {(\mathit{R},\mathit{V})}_{\mathrm{b}/\mathrm{t}},\\ {(\mathit{R},\mathit{V})}_{\mathrm{t}},\\ m_{\mathrm{p}},m_{\mathrm{i}},m_{\mathrm{o}},m_{\pi }\end{array}\right)& \leftrightarrow \hfill & \hfill \left(\begin{array}{c}{a}_{\mathrm{p}},e_{\mathrm{I}},{\omega }_{\mathrm{I}},t_{\mathrm{asc}\mathrm{p}},i_{\mathrm{I}},{\mathrm{\Omega }}_{\mathrm{p}},\\ a_{\mathrm{b}},e_{\mathrm{O}},{\omega }_{\mathrm{O}},t_{\mathrm{asc}\mathrm{b}},i_{\mathrm{O}},{\mathrm{\Omega }}_{\mathrm{b}},\\ a_{\mathrm{t}},e_{\mathrm{\Pi }},{\omega }_{\mathrm{\Pi }},t_{\mathrm{asc}\mathrm{t}},i_{\mathrm{\Pi }},{\mathrm{\Omega }}_{\mathrm{t}},\\ m_{\mathrm{i}}/m_{\mathrm{p}},P_{\mathrm{I}},P_{\mathrm{O}},P_{\mathrm{\Pi }}\end{array}\right),\end{array}$$(3)

where a_x, e_X, ω_X, t_ascx, i_X, Ω_x are the semi-major axis, orbital eccentricity, argument of periastron, time of passage at the ascending node³, inclination of the orbital plane with respect to the plane of the sky, and longitude of ascending node of bodies (or effectively a centre of mass) x ∈ {p,b,t} and binaries X∈{I,O,Π}. We note that we use lower-case letters when an orbital element refers to a body in particular, while we use an upper-case letter when the orbital element is associated with the binary as a whole independently of its components. On the first three lines, the correspondence between state vectors and orbital elements is performed using the post-Newtonian formalism of Damour & Deruelle (1985) (see also Voisin et al. 2020b). Except for the inner binary, this is equivalent to using the usual definition of orbital elements (e.g. Beutler 2004). On the last line, the masses m_{p, i, o, π} are related to the mass ratio m_i/m_p and orbital periods P_I, O, Π using Kepler’s third law. The state vectors of each individual bodies are derived from the left-hand side of Eq. (3) using centre-of-mass relations at first post-Newtonian order Voisin et al. (2020b) (the centre of mass of the whole system being at the origin of coordinates, by definition).

In practice, the fitted parameters are combinations of the orbital elements in Eq. (3). These choices reflect usual practices in pulsar timing, as well as relative sensitivity and correlations between parameters (see Tables E.1 and 1). Of particular interest are the use of Laplace-Lagrange parameters e_Xcosω_X, e_Xsinω_X particularly adapted to low-eccentricities (e.g. Lange et al. 2001) and the use of δi = i_I − i_O, δΩ = Ω_b − Ω_p, which allowed us to better reflect the fact that the two orbits (I, O) are co-planar within error bars (see below).

3.2.2. Effects of mutual interactions

Anticipating the results of Sect. 4, the only measurable effect of a small planet is through the induced Rømer delay. This results from the fact that its orbit is not relativistic, and its small mass and large orbital inclination suppress the Einstein or Shapiro delay down to the nanosecond level at most. The Rømer delay is caused by the variation of the projected distance of the pulsar onto the line of sight of an observer located at the Solar System barycentre,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_R=-\frac{{\mathit{n}}_{\odot }·{\mathit{R}}_{\mathrm{p}}}{c}=-\frac{1}{c}({\mathit{n}}_{\odot }·{\mathit{R}}_{\mathrm{p}/\mathrm{b}}+{\mathit{n}}_{\odot }·{\mathit{R}}_{\mathrm{b}/\mathrm{t}}+{\mathit{n}}_{\odot }·{\mathit{R}}_{\mathrm{t}}),\end{array}$$(4)

where n_⊙ is the unit vector pointing from the Solar System to the pulsar system barycentre, and the right-hand side is obtained using Eq. (2).

In the present analysis, R_p is computed numerically. However, it is interesting to use the decomposition of Eq. (2) in order to estimate what can be measured. We show in Appendix B that the three terms in Eq. (2) map to Jacobi coordinates. These coordinates permit a perturbative treatment of the orbital dynamics, which shows that, as said above, each term follows a Keplerian motion at leading order. That means that in first approximation, each term in Eq. (4) can be expressed using the usual binary expression of the Rømer delay (see Eq. (B.16) or e.g Hobbs et al. 2006; Lyne & Graham-Smith 2012). However, the sole detection of the Keplerian Rømer delay only allows one to constrain five parameters out of seven associated with the planetary binary: the orbital period P_X, the semi-major axis projected onto the line of sight a_xsini_X, the orbital eccentricity e_X and argument of periastron ω_X, where x and X are defined as above.

As per our parametrisation, the mass of the planet is the solution of Kepler’s third law:

$$\begin{array}{c}\hfill a_{\mathrm{t}}^3{n}_{\mathrm{\Pi }}^2=G\frac{m_{\pi }^3}{{(m_{\mathrm{t}}+m_{\pi })}^2},\end{array}$$(5)

where n_Π = 2π/P_Π is the planet’s mean motion. The three stellar masses, and therefore the mass m_t = m_p + m_i + m_o, can be determined thanks to mutual interactions, Shapiro and Einstein delays within the triple system (Voisin et al. 2020b). However, a_x is degenerate with sini_X at Keplerian order.

In Appendix B, we derive a simplified model of the first-order correction to the orbital motion of the planetary binary using Laplace-Lagrange perturbation theory. It amounts to decomposing ${\mathit{R}}_{\mathrm{t}}={\mathit{R}}_{\mathrm{t}}^{\mathrm{K}}+\delta {\mathit{R}}_{\mathrm{t}}$ where the first term describes a Keplerian orbit of the planet with a point mass m_t at t, and the second term accounts for the finite size of the triple system to first order in the ratio between the size of the planetary binary and that of the outer binary. It leads to a Rømer delay term δΔ_R = −n_⊙ ⋅ δR_t/c of the form, Eq. (B.22),

$$\begin{array}{c}\hfill \delta {\mathrm{\Delta }}_R=\alpha t({\gamma }_{c}cos{n}_{\mathrm{O}}t+{\gamma }_{s}sin{n}_{\mathrm{O}}t),\end{array}$$(6)

where n_O = 2π/P_O is the mean motion of the outer binary and α, γ_c, γ_s are constants defined in Eqs. (B.23)–(B.25). This functional time dependence, a sinusoidal oscillation at the outer-binary period the amplitude of which is growing linearly in time, is not degenerate with other components of the timing formula and can in principle allow for the separate measurement of αγ_c and αγ_s.

Importantly, γ_c and γ_s are independent functions of inclination i_Π and longitude of ascending node Ω_t, while α is derived from these as well as the Keplerian parameters of the system. This means that it is sufficient to detect the leading order correction of Eq. (6) to constrain all seven parameters associated with the planet’s orbit and mass.

3.3. Test of the equivalence principle

Each of the above two model categories is also divided in two depending on whether GR is assumed to be correct or if GWEP is being tested. The equivalence principle can be parametrised at Newtonian and post-Newtonian order by an effective, body-dependent, gravitational constant characterising the interaction between two bodies a and b, that is G_ab = G_ba = G(1+Δ_ab), where G is the Newtonian gravitational constant (see also Sect. 1 and Eq. (20) of Voisin et al. 2020b). In GR, Δ_ab = 0.

Solar-system experiments put strong constraints on the weak equivalence principle, as well as the strong equivalence principle in the weak-field regime (Touboul et al. 2019; Hofmann & Müller 2018). As shown in Voisin et al. (2020b), we can use these constraints to assume that within the precision of our experiment and (at least) within the framework of Bergmann-Wagonner scalar-tensor theories, any detectable violation of SEP can only be due to the pulsar, which is the only strongly self-gravitating object. Thus, the only sensitive parameter is Δ ≡ Δ_pb, where b ∈ {i, o, π}.

4. Results

We have compared seven different models all assuming that GR is correct. They are summarised in Table 2. We ran two additional models testing for SEP violation (see below). We inferred the posterior distribution functions (PDF) of each model using a Markov Chain Monte Carlo algorithm (MCMC) following the same procedure as in Voisin et al. (2020b). In particular, we used the same astrometric priors as in Voisin et al. (2020b) derived from Gaia DR2 (Lindegren et al. 2018) and spectroscopic observations of the inner white dwarf (Ransom et al. 2014; Kaplan et al. 2014). We used our own implementation (Voisin et al. 2020b) of the affine-invariant ensemble-sampling algorithm described in Goodman & Weare (2010); Foreman-Mackey et al. (2013).

The analysis was carried out on a chain sample with a length of at least 60 000, corresponding to more than 100 ensemble autocorrelation times. We used a set of 192 walkers saved to the chain every five iterations of the stretch move (see Goodman & Weare 2010). We evaluated the accuracy of derived statistics (e.g. Dunkley et al. 2005) by estimating the standard deviation of the mean and standard-deviation estimators, ${\widehat{\sigma }}_{\widehat{m}}$ and ${\widehat{\sigma }}_{\widehat{\sigma }}$ respectively, over at least 50 sub-chains (each spanning at least two ensemble autocorrelation times). The ratios of these quantities to the full-chain standard deviation σ, ${\widehat{\sigma }}_{\widehat{m}}/\sigma$ and ${\widehat{\sigma }}_{\widehat{\sigma }}/\sigma$ respectively, are used as ‘convergence ratios’ (Dunkley et al. 2005)⁴. We considered convergence was reached when both indicators were below 0.1 for every parameter chain.

4.1. Overall model comparison

We compared the best-fitting solutions of each model using both the Akaike (AIC) (Akaike 1974) and Bayesian Information Criteria (BIC) (Schwarz 1978) (see Table 2). To do so we computed the best-fitting solutions of each model by fitting them to the data using the deterministic minimizer Minuit (James & Roos 1975) starting from the best-fitting solution derived from the MCMC run. The obtained solution was always very close to the MCMC one, as expected, but important compared to the modest differences we found between models.

We ran four PLn models, from n = 2 to n = 5 Fourier components, initialising the PLn + 1 run from the outcome of the PLn run in order to limit computational cost. We ran two planetary models: a fully Keplerian one (called ‘Kepler’), and a numerical integration (‘Planet’)⁵ including all the effects of mutual interactions. We note that the PL2 model is equivalent to the Kepler model at leading order in eccentricity (see Appendix A), which is why we denote it ‘PL2/Kepl1’. Given the small difference between Kepler and PL2/Kepl1, no additional MCMC was done for Kepler, and its best-fitting solutions are based on a Minuit refit of PL2/Kepl1, thus saving on computing power.

Assuming the presence of a planet, the PL2/Kepl1 and Kepler models differ mostly by the addition of a third harmonic corresponding to the second-order term in eccentricity of the Keplerian Rømer delay (see Appendix A). Contrary to the PL3 model, the phase of the third harmonic is not independent from the two others, and its amplitude is not determined following a power law. Given e ∼ 0.25 and a projected semi-major axis x ∼ 4 μs (see Table 1), the amplitude of this harmonic is only about ○(xe²)∼0.25 μs. However small, this additional component appears to be captured by the Kepler model as both AIC and BIC are significantly better than for PL2/Kepl1 in Table 2. This can be visualised in Fig. E.1 where the small but significant peak (< 0.5% probability compared to white noise) at the third harmonic of the PL2/Kepl1 model residuals disappears with the Kepler model (see also Fig. 1). Thanks to its limited number of parameters the Kepler model has the best BIC of all.

Fig. 1. Comparison of numerically computed Rømer delay with first-order and Keplerian-order approximations. Approximate expressions have been least-square fitted to the numerical result. Top: Three versions as well as residuals of the least-square fit (lower panel). Bottom: Lomb-Scargle periodogram of the three versions. Vertical dashed lines mark the fundamental (black), second harmonic (grey), and third harmonic (dash-dotted grey) of the planet period. Vertical dotted lines mark the fundamental (black) and first harmonic (grey) of the outer-binary period.

The Planet model has a moderately better χ² but due to its two additional parameters accounting for inclination and longitude of ascending node, its AIC and BIC are not as good as those of the Kepler model, and even worse than PL2/Kepl1 concerning BIC. This suggests that mutual interactions may be detected, since the χ² is better, but only marginally, which is in line with the large uncertainty on the planet mass (see Table 1).

Assuming the red-noise hypothesis, that is the PLn models, PL3 is clearly favoured by BIC and only marginally worse than PL5 according to AIC (ΔAIC ≃ −1.38). Thus, we make it our reference achromatic red-noise model in the following. PL3 has the best χ² of all models, moderately better than Planet (Δχ² ≃ −2.4), and a moderately better AIC than Kepler (ΔAIC ≃ −2.0) but significantly worse BIC (ΔBIC ≃ +5.4).

In Table 2, one may notice that PL4 is actually a worse fit than both PL3 and PL5. In theory, since PL3 is nested in PL4, the χ² of the latter is expected to be at least as good as that of the former. We interpret the fact that it is not the case here as an occurrence of one of the local minima, which we routinely encountered while fitting. More details are given in Appendix C, but it is sufficient to say here that although Table 2 is the result of our best efforts within the limits of reasonable computing resources, it cannot be excluded that the differences between the models are below the level of systematics induced by local minima.

Fig. 2. Dispersion measure per time interval in the model PL3DM10. The x axis gives the time interval index, and the intervals are equal. Error bars delimit the 68% confidence region.

We also ran a model called PL3DM10, for which the MCMC only sampled the six achromatic red-noise parameters, the ten DM_X DM parameters as well as the pulsar spin parameters $\overline{f}$, $\overline{·}f$ (see Table E.1) since they may also correlate with the long-term effects of red noise. All orbital and astrometric parameters were fixed to their best-fitting PL3 values so as to optimise computational costs. We note that DM and DM′ were also fixed to their PL3 value, thus subtracting that linear trend from the computed DM_X. The resulting DM values are shown in Fig. 2 as a function of the corresponding time interval. No clear DM variation is seen given that the dispersion of the DM_X values is comparable to their uncertainties, and to the uncertainty on the global DM parameter as well (Table E.1). However, we note that AIC marginally favours PL3DM10 over PL3 (ΔAIC = −1.3) but that BIC strongly rejects it (ΔBIC = 58) owing to the large number of additional parameters required.

4.2. Details of the Planet and PL3 models

The inherently different nature of the Planet and PL3 models justify a more detailed look at both of them. If the planet hypothesis is restricted to the Kepler model, which can be seen as a sub-model of Planet, then it performs similarly to PL3 across all statistical metrics in Table 2. In addition, the potential systematic errors in best-fit estimations do not reasonably allow us to favour one hypothesis over another solely based on the statistical criterion (Appendix C).

We have collated the results of the MCMC Planet and PL3 runs in Table E.1 for parameters common to the two models and Table 1 for model specific parameters. Additionally partial correlation plots of each PDF are given in Appendix D, and complete ones can be found on the repository in footnote².

The timing signal introduced by both models are shown in Fig. 1. In that figure we have also represented the signal from the PL2/Kepl1 model, the Kepler model (Appendix A), as well as the perturbative model described in Appendix B. It can be seen that the Kepler model captures the third orbital-frequency harmonic similarly to PL3, while the perturbative model captures the dominant contribution from mutual interactions but lacks power at the third harmonic due to its limitation to first order in eccentricity. In the frequency domain, one sees that the largest difference between the PL3 and Planet models peaks at the outer binary frequency, as expected from the perturbative model, and decays in lower and upper frequency tails around that peak. However, as can be seen in the periodogram of best-fit residuals in Fig. E.1, this particular frequency does not allow for a clear signature with the current data as it is already very well fitted due to the outer binary contribution. In the case of PL3, Kepler, and PL2/Kepl1, it even appears to be over-fitted. Besides that, both the Planet, Kepler, and PL3 models appear consistent with white-noise in Fig. E.1 insofar as their residuals lie within the 95% confidence region for such noise.

In Table E.1, parameters that are common to both models are equal within error bars. Uncertainties on orbital parameters appear to be similar although systematically slightly smaller in the Planet model. On the other hand, uncertainties of spin parameters are an order of magnitude smaller in the PL3 model. These two parameters are the only two that significantly correlate with the low-frequency signal due to their ability to effectively absorb parabolic trends. Therefore, this difference probably comes from the two different ways that spin parameters correlate with model-specific parameters.

It is noteworthy that the orbital plane of the inner binary still cannot be separated from that of the outer binary, which is reflected by their relative inclination δi and relative longitude of ascending nodes δΩ being consistent with zero within error bars in both models. We also note that the tension in declination proper motion μ_δ compared to its Gaia prior that was observed in Voisin et al. (2020b) has disappeared in the present results.

The longitude of ascending node Ω_b is incompatible with that of Voisin et al. (2020b). During the course of the MCMC runs (irrespective of the models), a peak at the beat frequency between the Earth orbital period and outer-binary period became clearly visible in the residual periodogram as the low-frequency signal was being fitted out. This peak is characteristic of the annual-orbital parallax (Kopeikin 1995), which in the absence of mutual interaction is the only way to measure Ω_b by timing⁶. However, the MCMC seemed unable to find a solution accounting for this peak while staying in the vicinity of the initial value of Ω_b, and the assumption was made of a local minimum around the Ω_b value from Voisin et al. (2020b). The MCMC chain was restarted with values of Ω_b offset by ±90deg and 180deg from that of Voisin et al. (2020b). All converged to the value reported in Table E.1, which successfully fits out annual-orbital parallax. This value is 180deg away from Voisin et al. (2020b) (within uncertainty), thus suggesting that mutual interactions in the triple system can constrain the longitude of ascending node with a degeneracy of ±180deg lifted by the detection of annual-orbital parallax. Dedicated analytical or numerical work is needed to verify this conjecture.

An important point is that the uncertainty re-scaling parameter EFAC is now down to ≃1.11 from ≃1.31 in Voisin et al. (2020b), which translates into a significant improvement in accuracy. This is due to the fact that EFAC was previously absorbing the dispersion resulting from the unaccounted low-frequency signal. On the other hand, a fit with Minuit (James & Roos 1975) indicates that EFAC would need to be up to ∼1.6 in order to accommodate that signal in the longer dataset of this paper if it remained unmodelled

4.3. Galactic motion

For J0337 we are in the rare position of knowing the 3D velocity of a pulsar system with respect to the Solar System barycentre. The transverse velocity (magnitude and direction) is obtained from proper motion and distance (see Table E.1). High resolution spectroscopy of the inner white dwarf allowed for the determination of the systemic radial velocity: 29.7 ± 0.3 km s⁻¹ (Kaplan et al. 2014). With this information at hand, starting with the current location of the pulsar one can integrate its Galactic motion back in time. To do this, we use the Solar position and velocity parameters, the Galactic potential, and the software of McMillan (2017). Figure 3 shows the Galactic orbit of the J0337 system during the past 500 Myr. It is obvious that the system follows approximately the overall Galactic rotation. Indeed, its speed relative to the frame co-rotating with the Galaxy remains in the range ∼[30, 55] km s⁻¹ in the entire integration span (see Fig. 4). This finding is of particular importance for Sect. 5.3, where we discuss the evolutionary history of the system, as it suggests that the formation of the pulsar had only a small kick imparted on the system.

Fig. 3. Galactic motion of the J0337 system during the past 500 Myr. The blue dot marks its current position. The orange dot shows the location of the Sun. The orbit was calculated with the Galactic gravitational potential and the software provided by McMillan (2017).

Fig. 4. Velocity of the J0337 system with respect to the local galactic co-rotating frame during the last 500 Myr.

4.4. Tests of the strong equivalence principle with the Planet or PL3 models

Fig. 5. Measurements of the SEP violation parameter Δ (left column) and its absolute value |Δ| (right column) assuming the PL3 (upper row) or Planet (bottom row) models. Vertical dotted lines mark the mean value of the distributions, while the vertical dashed lines delimit the 95% confidence regions.

In order to test the SEP, we completed MCMC runs of the PL3 and Planet models while letting the SEP Δ parameter free (instead of fixed Δ = 0 when assuming GR). For other parameters results are compatible with the GR runs reported in Table E.1 and 1, although with much wider error bars for the orbital parameters due to their strong correlations with Δ. Besides the extended dataset and lower EFAC, this mainly explains why uncertainties reported in Table E.1 can be orders of magnitudes smaller than in Voisin et al. (2020b), which reported only an SEP run. In Fig. 5, we report the marginal PDF of the SEP Δ parameter for both models. Their 95% confidence region can be expressed as

$$\begin{array}{cccc}\hfill \mathrm{\Delta }& =0.{97}_{-1.48}^{+1.58}\times {10}^{-6}\hfill & \hfill \text{and}& |\mathrm{\Delta }|<2.29\times {10}^{-6}\text{(PL3)},\hfill \end{array}$$(7)

$$\begin{array}{cccc}\hfill \mathrm{\Delta }& =0.{51}_{-1.28}^{+1.14}\times {10}^{-6}\hfill & \hfill \text{and}& |\mathrm{\Delta }|<1.46\times {10}^{-6}\text{(Planet)}.\hfill \end{array}$$(8)

In both cases the width of the confidence region of Δ is reduced by ∼12% with respect to Voisin et al. (2020b). This can be related to the reduction by ∼16% of the EFAC parameter assuming a locally linear dependence on this parameter. However, the bound on |Δ| is improved by 30% with respect to Voisin et al. (2020b) in the Planet model while it worsens by 10% in the PL3 model. The latter case is due to the fact that the mean value of the distribution is offset from zero by more than one sigma.

5. Discussion

Our results show that DM variations cannot be responsible for the ∼4 μs amplitude of the observed low-frequency signal (Fig. 2), but may at most contribute marginally at the ∼0.1 μs level, if at all. Thus, in what follows we consider only the achromatic models.

The red-noise model following a power-law spectrum with three Fourier components (the PL3 model) provides the best description of the data from an information-theoretic point of view. However, the physical nature of the signal in the planet and the red-noise hypotheses are completely different, while the difference in terms of fitting residuals is minute, as can be seen in Figs. E.1 and 6. Besides, even in the Planet model a moderate red-noise component is not unexpected, typically at the ∼0.1 μs level. We did not try such a planet+red noise model since the additional complexity seems disproportionate with respect to the Planet model residuals, and would likely lead to over-fitting. However, our main goal was to address the large 4 μs signal, which both red-noise and Planet models do successfully. Thus, in what follows we discuss separately the two hypothesis.

Fig. 6. Timing residuals. Top: Best-fit residuals of model PL3. Bottom: Difference between residuals of best-fit PL3 and Planet models. Error bars are not shown for clarity, but the median 1-sigma uncertainty is 1.9 μs, and the mean is 2.2 μs.

5.1. Achromatic red noise hypothesis

Red noise in millisecond pulsars has been extensively studied in the frame of pulsar timing arrays (PTAs). Here we compare the red-noise properties of J0337 to those of the 67 different pulsars⁷ studied by the NANOGrav (Alam et al. 2021), EPTA (Chalumeau et al. 2022), and PPTA (Goncharov et al. 2021; Reardon et al. 2021). In those works, red noise is modelled as a Gaussian process over a Fourier basis characterised by a power spectrum density of the form P(ν)=αA_GP²(ν/yr⁻¹)^−γ_GP where α = 1/12π² in Chalumeau et al. (2022); Goncharov et al. (2021) or α = 1 in Alam et al. (2021), A_GP is the amplitude of the process at frequency ν = 1 yr⁻¹ and γ_GP is its exponent. Although we did not use a Gaussian process framework and therefore cannot make a rigorous comparison, we may estimate the parameters of a Gaussian red-noise process that would produce the same average spectrum as the one we fit in this paper. Up to an order one factor, we get A_GP ∼ A(T/α)^1/2(ν/yr⁻¹)^γ_GP/2 where we used the notations of Eq. (1), T is the time span of observations, and γ_GP ∼ 2γ. Using the results of our PL3 model, we obtain⁸ $A_{\mathrm{GP}}\sim {\alpha }^{-\frac{1}{2}}0.04\mathrm{\mu }\mathrm{s}{\mathrm{yr}}^{1/2}$ and γ_GP ∼ 5.4.

The value of the exponent places it among the steepest red-noise spectra, as is expected from a very smooth, quasi-sinusoidal signal. On the other hand, the amplitude of the signal, ∼4.4 μs at a period of 8 years (∼4.4 μs@8 yr), is only surpassed by PSR J1824-2452A with ∼8 μs@8 yr and γ_GP ≃ 5 among pulsars with a similarly steep spectrum (Goncharov et al. 2021). It is followed by PSR J1939+2134 (PSR B1937+21) with ∼1 μs@8 yr and γ_GP ≃ 5.4 (Goncharov et al. 2021; Reardon et al. 2021). Interestingly, PSR J1824–2452A and PSR J1939+2134 are both somehow extreme MSPs. They have the second and third largest spin-down power of all MSPs, $\dot{E}\sim 2\times {10}^{36}\mathrm{erg}/\mathrm{s}$ and ∼1 × 10³⁶ erg/s respectively (Reardon et al. 2021), as well as the second and fourth largest spin-down rates, according to the ATNF catalogue (Manchester et al. 2005), with $\dot{f}\sim -2\times {10}^{-13}\mathrm{Hz}/\mathrm{s}$ and $\dot{f}\sim -4\times {10}^{-14}\mathrm{Hz}/\mathrm{s}$ respectively Reardon et al. (2021). They are also among MSPs with the youngest characteristic ages with 30 and 240 million years, respectively. All these quantities are one to two orders of magnitude above the bulk of MSPs. Both pulsars are also among the few MSPs known to produce giant pulses (e.g. Knight et al. 2006). These characteristics are significant whether one considers that achromatic red noise is due to activity in the magnetosphere (Lyne et al. 2010) or to turbulence in the superfluid core of the neutron star (Melatos & Link 2014) since in both cases a correlation with spin-down rate is expected.

Empirically, a correlation is indeed observed over the general pulsar population (Hobbs et al. 2010; Lyne et al. 2010; Shannon & Cordes 2010). A scaling law $\propto f^{\alpha }{|\dot{f}|}^{\beta }$ can be fitted to the amplitude of red noise. In particular, Shannon & Cordes (2010) found α ≃ −1.4; β ≃ 1.1. At that time only a few MSPs showed detectable noise, notably PSR J1939+2134, and it appeared consistent with that scaling law. Restricted to a population of pulsars with comparable spin frequencies f this law means that noise should correlate with spin-down rate $\dot{f}$ or even spin-down power. More generally, it indicates that the characteristic spin-down age $\tau =f|\dot{f}{|}^{-1}$ can be used as a reasonable estimate, as was also noted in Hobbs et al. (2010): the younger the more noisy.

On the other hand, J0337 does not show any uncommon intrinsic properties but a spin-down power of $\dot{E}\sim 3\times {10}^{34}\mathrm{erg}/\mathrm{s}$, spin-down rate of $\dot{f}\sim -2\times {10}^{-15}\mathrm{Hz}/\mathrm{s}$, and characteristic age of 2.4 billion years. These quantities place J0337 in the bulk of the distribution of MSPs and according to the aforementioned scaling law its noise amplitude should be at least five to ten times weaker. We note that this is actually the case of the other MSPs with comparably steep index (> 4) in Chalumeau et al. (2022); Goncharov et al. (2021); Alam et al. (2021). However limited, these comparisons suggest that if indeed red noise causes the observed signal from J0337 then it is unusually strong, especially when one considers how steep its spectrum is.

In some of the occurrences where the timing noise is strong, a correlated variation of the pulse profile has been observed (Lyne et al. 2010). Consequently, we have tested the presence of such a slow variation of the pulse shape. To do so, we have built a series of 100-day integrated profiles distributed over the whole observation span. We then used the pat tool of PSRCHIVE in order to compare them to our main template used for ToA determination after correcting for phase shifts. The differences between the main template and sub-integrations can be quantified by calculating their reduced χ², as shown in Fig. 7. The reduced χ² is usually close to 1, indicating that we do not detect any significant pulse shape variations, except in May–Oct. 2019 for ∼150 days (MJD 58631 to MJD 58780). This corresponds to the period of instrumental modifications (see Sec. 2) that we have eventually excised. In order to quantify the dispersion of reduced χ², we show their two-standard-deviation interval in Fig. 7 computed on the clean χ² sample only. All but one element fall into that region, in line with expectations. Indeed, each profile possesses ∼2000 degrees of freedom, thus the χ² distribution is well approximated by a Gaussian law and the above-mentioned interval corresponds to a probability of 95%.

Fig. 7. Reduced χ2 of the residuals between each 100-day sub-profile and the main template as a function of date. Two sub-profiles are in the time interval when imperfect polarisation calibration was performed (‘UnperfectPolCal’, see main text). The inset shows the residuals obtained for one of these sub-profiles, while others do not show discernible structure (i.e. noise). The two horizontal lines delimit the two-standard-deviation interval around the mean of the reduced χ2 distribution corresponding to σred.χ2 ≃ 0.043 (imperfect polarisation region excluded).

5.2. Planet hypothesis

Assuming a mere Keplerian orbit of the planet without mutual interactions, the Kepler model appears to be favoured by BIC and moderately disfavoured by χ² and AIC. The full Planet model does not improve sufficiently its χ² in order to improve either AIC or BIC compared to Kepler, but its two additional parameters are nonetheless constrained by MCMC albeit with relatively large uncertainties. This allows us to discuss the complete solution constraining all seven parameters reported in Table 1, and in particular what they indicate about the formation and stability of the system.

The Planet model points to a very low-mass planet in a wide eccentric orbit around the stellar triple system. The mass of the planet is $1.{23}_{-0.66}^{+1.1}\times {10}^{-8}{M}_{\odot }$, which is about 1.3 × 10⁻⁵ M_Jup, ∼0.004 M_Earth or ∼0.4 M_Moon. As such, it is lighter than any of the planets orbiting PSR B1257+12 (Konacki & Wolszczan 2003) and could be the lightest exoplanet to date according to the extrasolar planet encyclopedia⁹.

We have evaluated the short-term stability of the orbit of the candidate planet by applying the frequency analysis method of Laskar (1988); Laskar (1990); Laskar (1993); Laskar (2005) in the eccentricity-period plane, Fig. 8, with the help of the TRIP software (Gastineau & Laskar 2011). These two parameters are relevant for stability, as they control the distance to the triple system, and other parameters are fixed to their best-fitting values. If the system was perfectly integrable, it would possess an intrinsic set of fundamental frequencies (related to the existence of action-angle coordinates), which would be constant over time. These frequencies could be retrieved very accurately by running the frequency analysis method on a numerical solution of the system. For a weakly chaotic system, however, fundamental frequencies can still be defined, but they are now valid only in a restricted interval of time, as chaos makes the system slowly diffuse in the available region of the phase space. The measure of the drift of the fundamental frequencies then gives a direct indication of the level of chaos present in the system (see Laskar 1990). Here, we are interested in the fundamental frequency n of the planet associated with its orbital motion (i.e. its ‘mean’, or ‘proper’ mean motion): its drift rate informs us about the short-term orbital stability of the planet. In practice, we run the frequency analysis on the two halves of a 2000-year numerical integration of the orbits; Fig. 8 shows the difference between the two values of n obtained.

Fig. 8. Chaos map representing the variation of the proper mean motion |nf − ni| between two halves of a 2000-year integration on a grid of 120 × 120 initial conditions in the period-eccentricity plane of the planet. White pixels represent orbits that are so chaotic that we could not run the frequency analysis. The eccentricity grid initially went all the way to 0.9, but we show here only the 102 lower values, as nearly all pixels in the last 18 lines are white. Vertical dotted lines show the positions of resonances 8:1 to 12:1 with the outer-binary period. The cross shows the mean parameters of the planet as reported in Table E.1, and the ellipses represent the one and two-sigma confidence regions, respectively.

We can see in Fig. 8 that larger eccentricities or smaller periods lead to unstable motion due to the stronger influence of the inner triple system. Mean-motion resonances with the outer binary are also a noticeable source of chaos. Due to its large uncertainty, the confidence region of the planet overlaps with the unstable stripes associated with the 9:1 to 11:1 resonances. However, most of the region remains compatible with short-term stability, and it contains niches in which virtually no chaos is detected (values < 10⁻⁵). Due to the very short period of the inner binary, studying the long-term (Gyr) stability of the system would require building a dedicated averaged model, which is out of the scope of the present work. We have nonetheless pushed numerical integration of the best-fitting solution to 100 000 years with reasonable accuracy (and to 100 Myrs without GR), and we did not see any long-term drift of the planet’s orbital elements.

However, long-term numerical integrations of the best-fitting solution reveal an intriguing behaviour for the planet’s orbit. First, we note that the best-fit inclination of the planet’s orbital plane with respect to the plane of the inner two binaries is $\delta i_{\mathrm{\Pi }}={119}_{-42}^{+16°}$ (Table E.1). This indicates that the motion of the planet is very inclined, and very likely retrograde, with respect to the orbits of the two inner binaries. More puzzling, this confidence interval is centred on the inclination of the exterior von Zeipel-Lidov-Kozai resonance¹⁰, whose location can be approximated by $cos(\delta i_{\mathrm{\Pi }})=-1/\sqrt{5}$, that is, δi_Π ≃ 116.6° (see Gallardo et al. 2012; Saillenfest et al. 2016). Contrary to the classic von Zeipel-Lidov-Kozai resonance for an outer perturber (see von Zeipel 1909; Lidov 1962; Kozai 1962), this resonance appears beyond quadrupole order, so its width is quite narrow in inclination. And yet, long-term numerical integrations with initial conditions chosen near the best-fit solution show that the planet can stably librate within the resonance. This peculiarity is directly apparent when using a system of coordinates in which the third axis is aligned with the total angular momentum. In this reference frame, the orbital inclination of the two inner binaries is close to zero and the inclination of the planet is close to δi_Π. As the planet is in the Kozai resonance, its eccentricity and inclination are affected by correlated oscillations, and its argument of periastron oscillates around π/2 (instead of circulation between 0 and 2π see Fig. 9). As the resonance only spans a few degrees in inclination, obtaining this configuration in a hierarchical system with nearly circular inner perturbers may seem hard to attribute to chance – even though the confidence region is wider than the resonance.

Fig. 9. Evolution of the eccentricity and inclination of the planet as a function of its argument of periastron. Orbital elements are measured with respect to the orbital plane of the inner two binaries. The small dots show the trace of a 20-Myr numerical integration of the planet, as given by an N-body integrator without GR. Initial conditions are chosen near the best-fit orbit.

Considering this unlikely coincidence, we may conjecture that the Kozai resonance stabilises the planet’s orbit, possibly because of the bounded secular evolution of the argument of periastron ω_π it implies. Indeed, an argument of periastron confined near π/2 means that when the planet is at pericentre, it is always away from the orbital plane of the two binaries. If this conjecture is correct, then the planet may be the only remnant of a larger population of small bodies that were on less stable orbits, or alternatively may have formed there from the gas expelled during the recycling of the pulsar as discussed below.

5.3. Planet formation and survival

The formation of the triple compact-object system J0337 is truly remarkable, and modelling its formation is a major challenge (Tauris & van den Heuvel 2014). The system has survived at least three stages of mass transfer (including one common envelope (CE) phase) and the supernova (SN) explosion that created the pulsar. On top of this, the system has managed to remain dynamically stable on a long-term (Gyr) timescale. Tauris & van den Heuvel (2014) found a possible solution for J0337 using a self-consistent and semi-analytical approach. Their solution is not unique due to degeneracy, and it should thus be considered a solution rather than the solution. An interesting question is if their formation model can be extended to accommodate the formation, and survival, of a planet?

In the following, we adapt the formation model of Tauris & van den Heuvel (2014) – see their Table 1 and Fig. 1 for an overview of the various stages involved (with the corresponding stage numbers, which we refer to below). We begin by noting that in stage 2, it has been suggested that planets may form from the condensation of vast amounts of material ejected in the equatorial region during a CE phase (Beuermann et al. 2010; Schleicher & Dreizler 2014; Bear & Soker 2014).

Although we cannot rule out the possibility of an outer planet forming from material ejected by the giant-star progenitors of the WDs (stages 6 and 8), such a formation process seems unlikely. It would require a rare and complex mechanism, facing many challenges not present in traditional planet formation models (Drążkowska et al. 2023). One major issue is that, even for a single red-giant star, enough dust and gas would need to accumulate in the circumstellar envelope for planetesimals (the small building blocks of planets) to form. Additionally, the gas and dust would need to cool and condense over a prolonged period before a planet could develop, and this cooling phase is critical since high temperatures would inhibit planet formation. In the formation model of Tauris & van den Heuvel (2014), the material would likely be ejected via isotropic re-emission during RLO with high pressure and temperature, meaning the material would likely be blown away before planet formation could occur. That said, it remains a puzzle how some WDs and pulsars have orbiting planets, and further research is required on this topic.

Assuming a planet did form during the CE ejection in stage 2, we can still adopt the original scenario and model parameters derived by Tauris & van den Heuvel (2014), given the very small planet mass of only M_planet = 1.23 × 10⁻⁸ M_⊙. At stage 4, just prior to the SN explosion, we can calculate the probability that this planet survived the kinematic effects of the SN.

Fig. 10. Probability of a planet surviving the SN explosion as a function of its pre-SN orbital period and the present-day systemic velocity of J0337. The red circle marks our default value with $v_{\mathrm{sys}}=44{\mathrm{km}\mathrm{s}}^{-1}$ and $P_{\mathrm{orb}}^{\mathrm{planet}}=1000\mathrm{days}$, yielding Pbound = 33% (see text).

Figure 10 displays the probability of a planet surviving the SN explosion as a function of its pre-SN orbital period and the observed present-day systemic velocity of J0337. This probability can be calculated directly from the equations in Hills (1983):

$$\begin{array}{c}\hfill P_{\mathrm{bound}}=\frac{1}{2}\{1+\left[\frac{1-2\mathrm{\Delta }M/M-{(v_{\mathrm{sys}}/v_{\mathrm{rel}})}^2}{2(v_{\mathrm{sys}}/v_{\mathrm{rel}})}\right]\},\end{array}$$(9)

where M is the total mass of the system, ΔM is the amount of material ejected in the SN, $v_{\mathrm{rel}}=\sqrt{GM/a}$ is the relative velocity between the triple system (“body 1”) and the planet (“body 2”), where a is the separation between the centre of mass of the triple system and the planet (assumed to be in a pre-SN circular orbit), and finally where v_sys is the present-day observed velocity of J0337 with respect to its local standard of rest. The calculation is made by considering the triple system and the planet as an in-effect two-body problem, which is a good approximation given the large orbital separation of the planet and the fast-moving SN ejecta. In this picture, we therefore consider v_sys as the “kick” imparted on the triple system (“body 1”). For a general discussion of dynamical effects of asymmetric SNe in hierarchical multiple star systems, see Pijloo et al. (2012), and references therein.

The parameter values needed for Eq. (9) are as follows at stage 4: M = 4.10 M_⊙ (1.70 + 1.10 + 1.30 M_⊙), ΔM = 0.42 M_⊙, $v_{\mathrm{sys}}=44{\mathrm{km}\mathrm{s}}^{-1}$ and $P_{\mathrm{orb}}^{\mathrm{planet}}=1000\mathrm{days}$. The latter two values are our default values: $v_{\mathrm{sys}}=44{\mathrm{km}\mathrm{s}}^{-1}$ corresponds to our estimated value of the peculiar, that is, with respect to the local standard of rest, velocity of J0337 (see Sect. 4.3) and $P_{\mathrm{orb}}^{\mathrm{planet}}=1000\mathrm{days}$ is selected as a typical value from our range of solutions to the SN event based on Monte Carlo simulations with 2 × 10⁶ trials, following the recipe of Tauris et al. (2017). In these Monte Carlo simulations, we have imposed a criterion of a post-SN orbital period of the planet of 910 days (±3%). This value is found by calculating backwards from stage 9 (present day) to stage 5 (right after the SN explosion, and before the formation of the two WDs), considering the mass lost from the system when the main-sequence star progenitors of the WDs evolve to fill their Roche lobes and undergo mass transfer (Tauris & van den Heuvel 2023). The orbital separation ratio (before and after the mass loss) will simply scale inversely to the total mass ratio: a/a₀ = M₀/M. Based on our Monte Carlo simulations, we find maximum pre-SN orbital periods of {2430; 1890; 910; 410} days for v_sys={20; 44; 70; 100} km s⁻¹, respectively. The minimum pre-SN orbital period is about 100 days. It is limited by the condition of long-term dynamical stability of the system.

In terms of the effect of the interaction of the SN ejecta on the planet, the effect is expected to be negligible (Wheeler et al. 1975; Liu et al. 2015). This is mainly due to the large orbital separation of the planet at the moment of the SN. In addition, the surface area of the planet is presumably quite small, and its mass density is likely to be larger than that of a main-sequence star (∼1 g cm⁻³) or a gaseous planet (∼0.1 g cm⁻³). For comparison the Moon has a mean mass density of ∼3.3 g cm⁻³.

The large orbital inclination of $\delta i_{\mathrm{\Pi }}={119}_{-42}^{+16°}$ of the planet could have an origin in the SN explosion. It is certain that if the kick (or resulting recoil of the inner system) has a velocity vector that is not exactly in the orbital plane, then the system will be tilted. Whether or not this is sufficient to explain the orbital inclination of the planet is uncertain, and this inclination could also be due to dynamical interactions within the system.

To summarise, we find from calculations in Fig. 10 a most likely expected probability of the planet to survive the SN explosion of order 20%–50% (33% for our default value). Therefore, if a planet is present in the J0337 system with an orbital period of about 3310 days, and if this planet formed from ejected CE material in stage 2, then we find it quite likely that the planet could have survived the kinematic impact of the SN that created the neutron star (NS) in this system.

5.3.1. Constraining the NS kick

To place constraints on the kick velocity imparted directly on the NS at its formation is also possible, albeit rather complicated, and a proper treatment is beyond the scope of this paper. Using a step-by-step Monte-Carlo approach, however, where the kinematic effects are evaluated as in-effect two-body problems (step 1 consisting of “body 1”: NS+WD_i+WD_o and “body 2”: planet; step 2 consisting of “body 1”: NS+WD_i and “body 2”: WD_o; step 3 consisting of “body 1”: NS and “body 2”: WD_i) we find a rough constraint on the NS kick of w ≃ 110 − 125 km s⁻¹. Here we made the assumption that the direction of recoil (or “kick”) velocities is always isotropic. This is not likely to be the case because of the momentum in the pre-SN orbital plane in each step (only for the actual kick imparted onto the NS in step 3 is this a good assumption). Nevertheless, it seems to us that the resulting value of w may not change by much taking such effects into account. A NS kick value of w ≃ 110 − 125 km s⁻¹ is somewhat on the low side for measured pulsar kicks (with average values above 400 km s⁻¹, Hobbs et al. 2005). However, it is not an unusual small kick. Moreover, it is an obvious selection effect that the J0337 system exists, and thus the kick could not have been very large (see also discussions in Tauris & van den Heuvel 2014).

In addition the small kick inferred for J0337, and its relatively small velocity relative to its local galactic co-rotating frame (see Fig. 3) are very similar to those observed for PSR J1903+0327 (Freire et al. 2011, see especially their Fig. 10), a system that is thought to have originated as a triple, which later became unstable because of the widening of the inner binary (then a low-mass X-ray binary) caused by the mass transfer to the NS Freire et al. (2011); Portegies Zwart et al. (2011). This, again, emphasizes the importance of small SN kicks for the preservation of hierarchical systems

6. Conclusions

We have extended the work of Voisin et al. (2020b) by using a 700-day longer dataset of the Nançay timing data of PSR J0337+1715 and, most importantly, by attempting to model the pulsar’s residual long-term signal, which has an amplitude of ∼4.4 μs over ∼3000 days. To this end, we have complemented the numerical timing model for the triple stellar system introduced in Voisin et al. (2020b) with either an achromatic red-noise component or a planet in a hierarchical orbit around the triple system. Both models have been implemented in the code Nutimo (Voisin et al. 2020b). We note that the NRT ToAs considered in this study were extracted from total intensity profiles. In future work, we will consider using ToAs extracted from polarimetric profiles with the Matrix Template Matching (van Straten 2006) technique, which was shown by Guillemot et al. (2023) to significantly improve the quality of the NRT timing by compensating for imperfect polarisation calibration. (See footnote² for software and data release.)

We have shown in Sect. 4.2 that with sufficient data, the Planet model can display a clear signature that distinguishes it from red noise. Though we have attempted to select the best model using the AIC and BIC statistical criteria, we acknowledge that the difference between the models is within the systematic uncertainties in our numerical optimisation procedure given our current dataset. Nonetheless, both the red-noise and planet hypothesis are able to produce residuals that appear consistent with white noise (Fig. E.1).

The best red-noise model contains three Fourier harmonics, the amplitude of which decreases with a steep power-law index of ∼2.7 (or ∼5.4 in the power convention, see Sect. 5.1). The Planet model implies a Moon mass object in a hierarchical ∼3000-day orbit around the triple system, a ∼0.25 eccentricity, and a ∼119° inclination with respect to the fundamental plane of the triple system. Triple system, DM, and astrometric parameters are compatible within uncertainties such that they are essentially model independent. The full characterisations of the system assuming both models is presented in Tables E.1 and 1.

From a Bayesian point of view, one may wonder which model is a priori more probable. Planets around pulsars are rare (Kerr et al. 2015; Behrens et al. 2020; Nitu et al. 2022); however PSR J0337+1715 is itself a unique system with an unusual formation channel. Thus, it is possible that the probability of existence of a planet in such a system is not as small as in the general population. In Sect. 5.3, we have proposed a scenario derived from Tauris & van den Heuvel (2014) in which the planet forms from the material expelled during the first common envelope phase in the evolution of the triple system and survives the subsequent stages of evolution. In Sect. 5.2, we have shown the likely short-term stability of the inferred orbit and noted that the confidence interval on its inclination with respect to the triple system is centred on the value of the exterior von Zeipel-Lidov-Kozai resonance. This resonance is narrow, and the confidence interval still allows for a mismatch. However, if such a correspondence were confirmed, it seems unlikely that such a coincidence would occur. This leads us to conjecture that this particular orbit may be more stable, possibly making this planet the only remnant of a larger population formed out of the expelled material. As a by-product of this study, we estimate a rather low supernova kick of ∼110 − 125 km/s, not unlike the kick of PSR J1903+0327, which is believed to have originated as a triple system as well (Freire et al. 2011). This supports the idea that low kicks might be a condition of survival of such systems.

On the other hand, the achromatic red noise is a common phenomenon usually associated to processes intrinsic to the pulsar that are yet to be fully understood (Lyne et al. 2010; Melatos & Link 2014). However, the amplitude of the signal has been shown empirically to scale with quantities such as the spin-down age, spin-down rate, or power of the pulsar (Hobbs et al. 2010; Lyne et al. 2010; Shannon & Cordes 2010). In Sect. 5.1, we have shown that, according to those scalings, the red-noise amplitude necessary to explain the observed signal is a clear outlier by a factor of five to ten. Indeed, PSR J0337+1715 lies in the bulk of the millisecond pulsar distribution where red noise should be barely detectable with our data, if at all. A possible signature of strong red noise due to magnetospheric activity is its correlation to pulse profile variations, as observed in other pulsars (Lyne et al. 2010), but we did not find any evidence of such variations here. Contrary to the planet hypothesis, we currently do not see any reason to speculate that such a discrepancy could be related to the particular formation channel of the system.

We have updated our previous limit on violations of the strong equivalence principle (Voisin et al. 2020b). In Sect. 4.4, we have shown that the limit is somewhat model dependent, with |Δ|< 2.29 × 10⁻⁶ assuming the best red-noise model and |Δ|< 1.46 × 10⁻⁶ assuming the Planet model (both delimiting the 95% confidence region). If the latter represents a 30% improvement compared to Voisin et al. (2020b), the former is 10% worse. This systematic uncertainty emphasises the importance of selecting the correct model with future data.

Acknowledgments

This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This research has made use of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia. This work made use of the Scipy libraries (www.scipy.org). The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from the “Programme National Gravitation, Références, Astronomie, Métrologie (PNGRAM) funded by CNRS/INSU and CNES, France. We would like to thank the referee, Dr. Scott Ransom, for his insightful review, which allowed us to further improve this work.

References

Appendix A: Second-order Keplerian Rœmer delay

The Rœmer delay due to a pulsar orbiting in a binary system is given by (Blandford & Teukolsky 1976),

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{R}}=x(sin\omega (cosE-e)+\sqrt{1-e^2}cos\omega sinE),\end{array}$$(A.1)

where x is the projected semi-major axis , e is the orbital eccentricity, ω the argument of periastron, and E the eccentric anomaly. The latter is solution of Kepler’s equation,

$$\begin{array}{c}\hfill E-esinE=\sigma ,\end{array}$$(A.2)

where σ = n(t − T_p) with time t, time of passage at periastron T_p and n = 2π/P for orbital period P.

Kepler’s equation can be solved iteratively at any order in e using a fixed-point method following E_n + 1 = σ + esinE_n. Injecting E₂ into Eq. A.1 and Taylor expanding, we obtained the second-order Rœmer delay,

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{R}}& =x[sin\varphi \hfill \\ \hfill & +\frac{e}{2}(cos\omega sin2\varphi -sin\omega cos2\varphi -3sin\omega )\hfill \\ \hfill & -\frac{e^2}{8}\{(cos2\omega +4)sin\varphi -sin2\omega cos\varphi \hfill \\ \hfill & -3cos2\omega sin3\varphi +3sin2\omega cos3\varphi \}],\hfill \end{array}$$(A.3)

where we define ϕ = n(t − T_asc) with T_asc0 = T_p − ω/n, which corresponds to the time of ascending node at order e⁰. The first two lines correspond to the small-eccentricity timing model of Lange et al. (2001) where one can identify the Laplace-Lagrange parameters esinω and ecosω. In all rigour we left here the term in ∼ecosω, which is dismissed as an unnecessary additional constant in Lange et al. (2001) but can play a role in case of periastron precession (Susobhanan et al. 2018; Voisin et al. 2020a). Apart from this term, Eq. (A.3) was also derived in Zhu et al. (2019).

We can see in Eq. (A.3) that the effect of taking into account second-order terms is i) to modify the coefficients in front of first-harmonic terms and ii) to add signal at the third harmonic frequency. Without the third-harmonic contribution, i) would merely amount to an effective eccentricity and argument of periastron. However thanks to third harmonics this degeneracy is lifted and first-harmonic amplitudes become second-degree polynomials in eccentricity.

Most importantly in the frame of this work, we can see that once limited to first order, this timing model is equivalent to the PL2 model. Moreover, we see that all five elements x, e, ω, T_p, n can be measured independently at this order up to an error ○(e²). The equivalence between PL2 as given by Eq. (1), and this model approximately goes as follows:

$$\begin{array}{cc}\hfill n& =2\pi \nu ,\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill x& =A,\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill T_{\mathrm{p}}& =T_{\mathrm{ref}}-\frac{{\varphi }_2-{\varphi }_1}{\nu },\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill \omega & =2{\varphi }_1-{\varphi }_2,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill e& =2^{1-\gamma }.\hfill \end{array}$$(A.8)

The approximation lies in the fact that we neglected the difference between emission time and arrival time, Δ = t_a − t_e Einstein delay excluded, such that Eqs. (A.4)-(A.8) are valid up to corrections of order ○(Δ/P) <  < 1.

All subsequent harmonics are fully determined by these same five parameters, and therefore the Keplerian model has a higher predictive power than models such as PL3, which requires additional phase parameters. In addition it can be seen that contrary to the PLn models, amplitudes are not determined by a power law past PL2.

Appendix B: Simplified perturbative model

In order to assess how well the orbital motion of a small planet in a hierarchical orbit can be characterised we derive in this section an analytical model for a simplified configuration. The aim is to compute the leading order perturbation due to the planet beyond the Keplerian motion it induces on the barycentre of the triple stellar system (see below).

We model the system by three bodies: the first one corresponds to the inner binary, the second one to the outer white dwarf and the last one to the planet. The extent of the inner binary is thus neglected. We also assume Newtonian dynamics. We note that this three-body modelling could also be applied to the triple system itself since it is also hierarchical.

We denote M_0, 1, 2 the masses and ${\overrightarrow{R}}_{0,1,2},{\overrightarrow{V}}_{0,1,2}$ the position and velocity vectors. Indexes refer to the three bodies in the same order as above. The inner binary is represented by the position and velocity of its barycentre (R₀, V₀) and the sum of the masses of its two components. The justification for neglecting the dynamics of the inner binary is that due to its small extent relative to the separation with the planet it is much less sensitive to the planet perturbation than the outer binary. The problem is treated using Jacobi coordinates (e.g. Murray & Dermott 1999),

$$\begin{array}{cc}\hfill {\mathit{r}}_0& ={\mathit{B}}_2,\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill {\mathit{r}}_i& ={\mathit{R}}_i-{\mathit{B}}_{i-1}(i>0),\hfill \end{array}$$(B.2)

where the centre of mass of the first i + 1 bodies is

$$\begin{array}{c}\hfill {\mathit{B}}_i=\frac{1}{{\sigma }_i}\sum _{k=0}^i{M}_k{\mathit{R}}_k,\end{array}$$(B.3)

with ${\sigma }_i=\sum _{k=0}^i{M}_k$.

With these coordinates r₁ is the vector connecting the inner binary to the outer white dwarf, which we call the outer binary, and r₂ the vector going from the centre of mass of the triple stellar system to the planet. Given that the planet follows an orbit about five times wider than the size of the triple system and that its mass is negligible compared to those of the stars, we introduce the following scalings:

$$\begin{array}{ccc}\hfill ϵ& \sim \hfill & \hfill M_2/M_{0,1}\ll 1,\end{array}$$(B.4)

$$\begin{array}{ccc}\hfill \eta & \sim \hfill & \hfill r_2/r_1\sim R_2/R_{0,1}\ll 1,\end{array}$$(B.5)

where we use the notation $V=\Vert \overrightarrow{V}\Vert$.

Using canonical coordinates and their conjugate momenta P_0, 1, 2, the Hamiltonian of the system is given by

$$\begin{array}{c}\hfill H=\frac{1}{2}\sum _{k=0}^2\frac{{\mathit{P}}_k^2}{M_k}-\sum _{0\le k<l\le 2}\frac{G{M}_k{M}_l}{\Vert {\mathit{R}}_l-{\mathit{R}}_k\Vert },\end{array}$$(B.6)

where G is the gravitational constant.

We now decompose the Hamiltonian in an dominant integrable part H₀ and a perturbing term H₁ using Jacobi coordinates Eqs.(B.1)-(B.2) and their conjugate momenta p_i. The kinetic part of the Hamiltonian Eq.(B.6) keeps the same form after the substitution P_k → p_k, M_k → m_k where m₀ = σ₂ (the total mass of the system) and m_{i >  0} = m_iσ_i − 1/σ_i. Noting that p₀ is the total momentum of the system, we choose it to be zero without loss of generality. Further, the interaction term between zero and one is also left unchanged since r₁ = R₁ − R₀. On the other hand, the two terms involving the planet are expanded in powers of η. Collecting the terms,

$$\begin{array}{cc}\hfill H_0& =\frac{{\mathit{p}}_1^2}{2m_1}+\frac{{\mathit{p}}_2^2}{2m_2}-\frac{G{m}_1{\sigma }_1}{r_1}-\frac{G{m}_2{\sigma }_2}{r_2},\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill H_1& =-\frac{1}{2}\frac{{\sigma }_2}{{\sigma }_1}\frac{G{m}_1{m}_2}{r_2}[\frac{3}{2}{\left(\frac{{\mathit{r}}_1·{\mathit{r}}_2}{r_2^2}\right)}^2-{\left(\frac{r_1}{r_2}\right)}^2]+◯\left({\eta }^3\right).\hfill \end{array}$$(B.8)

One can see that H₁ is of order η² compared to the last term of Eq. (B.7) and of order ηϵ compared to the penultimate term, which justifies it being treated as a perturbation. We also note that σ₂/σ₁ = 1 + ○(ϵ), and therefore it is omitted in the following (however it would have to be kept if this methodology was applied to the triple stellar system).

The dominant Hamiltonian H₀ is recognised as the Hamiltonian of a mass m₁ orbiting around a central mass σ₁ and a mass m₂ independently orbiting around a central mass σ₂. Thus the solution is exactly given by two Keplerian orbits.

B.1. Timing model corresponding to the dominant Keplerian Hamiltonian

We consider only the largest delay, that is the Rœmer delay,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_R=-\frac{{\mathit{n}}_{\odot }·{\mathit{R}}_0}{c},\end{array}$$(B.9)

where n_⊙ is a unit vector along the direction going from the Solar System barycentre to the pulsar system barycentre and c is the speed of light. We note that the motion of the pulsar is approximated here again to the motion of the barycentre of the inner binary R₀. A complete treatment could formally be obtained by replacing R₀ by R_p = R₀ + δR_p where the second term accounts for the inner binary motion relative to its barycentre as well as every effect due to mutual interactions not accounted for in our treatment.

In terms of Jacobi coordinates,

$$\begin{array}{c}\hfill {\mathit{R}}_0=-\frac{M_1}{{\sigma }_1}{\mathit{r}}_1-\frac{M_2}{{\sigma }_2}{\mathit{r}}_2,\end{array}$$(B.10)

where according to the exact Hamiltonian H₀r_1, 2 follow Keplerian orbits of the form

$$\begin{array}{c}\hfill {\mathit{r}}_i=a_i{\left(\begin{array}{c}cos{E}_i-e_i\\ \sqrt{1-e_i^2}sin{E}_i\\ 0\end{array}\right)}_{F_i},\end{array}$$(B.11)

where F_i is the frame adapted to the orbital motion of coordinate i, that is with the direction of periastron along the x axis and orbital angular momentum along the z axis, a_i is the semi-major axis, E_i is the eccentric anomaly and e_i the eccentricity. Kepler’s equation relates E_i to the mean anomaly ℳ_i = 2π(t − T_i)/P_i, where P_i is the orbital period and T_i the time of periastron passage,

$$\begin{array}{c}\hfill E_i-e_{i}sin{E}_i={\mathcal{M}}_i.\end{array}$$(B.12)

In order to express the Keplerian Rœmer delay one needs to express both n_⊙ and R₀ in the same frame. We choose a frame where n_⊙ is the third basis vector. We call it the observer’s frame F_O (who would be standing at the Solar System barycentre). Orbital inclinations i_j are then defined relative to n_⊙ and longitudes of ascending nodes Ω_j define a rotation around n_⊙. Together with the arguments of periastron ω_j these angles define the rotations relating the observer’s frame and F_j,

$$\begin{array}{c}\hfill {\mathit{r}}_j{{|}_{F_O}=R_{{\mathrm{\Omega }}_j}{R}_{i_j}{R}_{{\omega }_j}{\mathit{r}}_j|}_{F_j},\end{array}$$(B.13)

where R_Ωj, R_ij, R_ωj are rotation matrices (see e.g. Beutler (2004)), and r_j|_F refers to the coordinate of r_j in frame F. From there, the formula expressing the Keplerian Rœmer delay Δ_j = −n_⊙ ⋅ r_j/c is well known (e.g. Lyne & Graham-Smith 2012; Hobbs et al. 2006),

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_j& =-\frac{a_{j}sin{i}_j}{c}\hfill \\ \hfill & (\sqrt{1-e_j^2}sin{E}_{j}cos{\omega }_j+cos{E}_{j}sin{\omega }_j-e_{j}sin{\omega }_j),\hfill \end{array}$$(B.14)

We note that the last term of Eq. (B.14) is usually omitted because it contributes only a constant delay to the timing solution, but we write it here for completeness as it can become important when orbits are perturbed (e.g. Voisin et al. 2020a). Let us also remark that this formula is independent of Ω_j and that a_j cannot be separated from sini_j, which leads to the impossibility of independently measuring these three parameters with the Keplerian Rœmer delay alone. Inserting Eq.(B.10) into (B.9), we can write

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_R=-\sum _{j=1}^2\frac{M_j}{{\sigma }_j}{\mathrm{\Delta }}_j.\end{array}$$(B.15)

B.2. Perturbation

We now turn to the perturbative treatment of H₁, Eq.(B.8), using Laplace-Lagrange’s variation of the orbital elements. In order to retain only the strongest effects and to keep the treatment as simple as possible we neglect all eccentricities, and we retain only secular terms. We note that we neglect eccentricity only in perturbative terms, not in leading order terms, that is Eq.(B.14). Thus we seek a Rœmer delay of the form Eq.(B.15) with

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_j={\mathrm{\Delta }}_j^{(0)}+\delta {\mathrm{\Delta }}_j,\end{array}$$(B.16)

where Δ_j⁽⁰⁾ is given by Eq.(B.14), and

$$\begin{array}{c}\hfill \delta {\mathrm{\Delta }}_j=\delta (-\frac{a_{j}sin{i}_j}{c}(sin(n_{j}t+{\varphi }_j))),\end{array}$$(B.17)

which denotes the variation of the leading order term in eccentricity of Eq.(B.14). For convenience we introduced ϕ_j = ω_j − n_jT_j such that at leading order n_jt + ϕ_j = E_j + ω_j.

The outer-binary term in Eq.(B.15) is dominant since r₂M₂/σ₂ ∼ (ϵ/η)r₁M₁/σ₁ and in the present case ϵ ≪ η. Therefore we focus on the perturbation of Δ₁. Without loss of generality we choose a frame such that Ω₂ = 0 such that Ω₁ is defined relative to the longitude of ascending node of the planet.

In order to get rid of periodic terms we average H₁ over ℳ₁, being understood that P₁ ≫ P₂, and obtain ${\overline{H}}_1={(2\pi )}^{-1}{\int }_0^{2\pi }{H}_1\mathrm{d}{\mathcal{M}}_1$. Although straightforward, this operation is very lengthy. We used the formal calculus software Sage (The SageMath Developers 2022) in order to manipulate the large expressions generated, without any particular difficulty. Therefore in what follows, we only report the important steps and the main results.

The resulting perturbing Hamiltonian does not depend on ω_1, 2, T₁, E₁, e_1, 2, Ω₂. The relevant Laplace-Lagrange variational equations are

$$\begin{array}{cc}\hfill {\dot{i}}_1& =-\frac{1}{n_1{a}_1^{2}sin{i}_1}\frac{\partial H_1}{\partial {\mathrm{\Omega }}_1},\hfill \end{array}$$(B.18)

$$\begin{array}{cc}\hfill {\dot{\mathrm{\Omega }}}_1& =\frac{1}{n_1{a}_1^{2}sin{i}_1}\frac{\partial H_1}{\partial i_1},\hfill \end{array}$$(B.19)

$$\begin{array}{cc}\hfill {\dot{\varphi }}_1& =-\frac{1}{n_1{a}_1}(\frac{1}{a_{1}tan{i}_1}\frac{\partial H_1}{\partial i_1}+2\frac{\partial H_1}{\partial a_1}),\hfill \end{array}$$(B.20)

$$\begin{array}{cc}\hfill {\dot{a}}_1& ={\dot{n}}_1={\dot{e}}_1=0.\hfill \end{array}$$(B.21)

We seek solutions of the form x(t)=x⁽⁰⁾ + δx(t) where x⁽⁰⁾ ∈ {a₁, n₁, i₁, Ω₁, ϕ₁} are the constant orbital elements that are solution of the unperturbed system. To leading order in perturbation we may thus substitute x → x⁽⁰⁾ in the right-hand side of Eqs.(B.18)-(B.21) and integrate in order to obtain $\delta x=\int \dot{x}\mathrm{d}t$. In order to focus on the main contributions we keep only terms diverging with time, that is terms ∝E₂. In the following, we drop the (0) exponent for readability since orbital elements always refer to the unperturbed solution, that is x = x⁽⁰⁾ unless otherwise stated.

Substituting x → x + δx(t) into Eq. (B.17) we obtained an expression of the form

$$\begin{array}{c}\hfill \delta {\mathrm{\Delta }}_1=\alpha t({\gamma }_{c}cos{E}_1+{\gamma }_{s}sin{E}_1),\end{array}$$(B.22)

where E₁ = n₁(t − T₁) at leading order in eccentricity and

$$\begin{array}{cc}\hfill \alpha & =\frac{3}{8}\frac{G{m}_2{a}_1}{c{n}_1{a}_2^3},\hfill \end{array}$$(B.23)

$$\begin{array}{cc}\hfill {\gamma }_s& =-sin{\mathrm{\Omega }}_{1}cos{i}_{1}sin{i}_2\hfill \\ \hfill & (cos{i}_{1}cos{i}_2+cos{\mathrm{\Omega }}_{1}sin{i}_{1}sin{i}_2),\hfill \end{array}$$(B.24)

$$\begin{array}{cc}\hfill {\gamma }_c& =sin{i}_1({cos}^2{i}_1{cos}^2{i}_2-\frac{10}{3})\hfill \\ \hfill & +cos{\mathrm{\Omega }}_{1}sin{i}_{2}cos{i}_{2}cos{i}_1(1+2{sin}^2{i}_1)\hfill \\ \hfill & +{(cos{\mathrm{\Omega }}_{1}sin{i}_2)}^{2}sin{i}_1(1+{sin}^2{i}_1).\hfill \end{array}$$(B.25)

In order to validate the ability of our perturbative approach to capture the main features of the Rœmer delay induced by the planet we have fitted Δ₁, obtained by substituting Eq.(B.22) in Eq. (B.16), to the numerically computed Rœmer signal. The fitting was done using a Levenberg-Marquardt method¹¹. The result is shown in Fig.1, where the Keplerian component Δ₁⁽⁰⁾ is also reported, and we can see that the perturbative formula successfully captures the dominating frequencies, that is the planet and outer binary orbital frequencies. The main residual appears at the first harmonic of the planet orbital frequency.

The parameters αγ_c, αγ_s can in principle be fitted independently due to the Fourier decomposition theorem. Since α only occurs as a common scaling in Eq.(B.22), we focus on γ_c, γ_s. If their relation to Ω₁, i₂ is not degenerate, Eqs. (B.24)-(B.25), then this allows the degeneracies of the Keplerian component to be lifted and full characterisation of the planet’s orbit and mass. In order to check for the absence of degeneracy, we plotted in Fig.B.1γ_c(Ω₁, i₂),γ_s(Ω₁, i₂) for i₁ = 39deg, which approximately corresponds to the observed inclination of the outer binary. We see that although not degenerate, γ_c is about five times larger than γ_s for this particular inclination i₁. Given the overall weakness of the signal, it might therefore be difficult to constrain accurately both i₁ and ω₂ in this particular case.

Fig. B.1. Maps representing the values of γc (left) and γs (right) as a function of (Ω1, i2) for i1 = 39deg.

Appendix C: Best-fit systematics

Table C.1 summarises all the fits that have been performed, while Table 2 only contains the best fit for each model. The procedure involved using either the mean, ‘mean, ’or the maximum-posterior probability (including prior), ‘max’ element of the MCMC sample of each model as the starting point for a deterministic fit using the Minuit library. Then, this ‘refit’ converges to the local optimum. The fact that the local optima generally depend on their starting points is evidence for a rough likelihood surface with multiple nearby local maxima. We note that each refit has been iterated one more time in order to make sure that this difference was not related to an approximate convergence of the algorithm.

Surprisingly, the PL4 model produced less optimal log-likelihood than the PL3 model, although the latter is nested in the former. In order to attempt to find a better solution we ran two MCMC: one starting from the last chain elements of the PL3 run and another from the PL5 run, which we call PL4bis. However, none of these fits performed better than PL3, suggesting additional complications in the probability landscape.

The Planet model behaves somewhat differently than the others. Indeed it is the one where the likelihood difference between the mean and max solutions is by far the largest, and the best fit is given by the max, suggesting that the posterior probability distribution is asymmetric or has multiple modes. This is consistent with its more irregular correlation plots, Fig. D.2, compared to the PL3 correlations in Fig. D.1. We also note that although the Planet max solution has the best log-likelihood of all models, the gain after refitting is more moderate than in other models.

To sum up, Table C.1 shows that a difference of a few units in log-likelihood may not be safely considered as significant given the systematic uncertainty of the fitting procedure.

Appendix D: Posterior distribution functions of model-specific parameters

Model specific parameters, as reported in Table 1, are not significantly correlated with other parameters of the timing model except for the re-scaled spin parameters $\overline{f}$, $\overline{·}f$. Therefore we show in this appendix the corresponding ‘correlation plots’ (or ‘corner plots’) showing the posterior distribution function of the PL5 and Planet model marginalised along all but two dimensions (or one across the diagonal), Fig. D.1 and Fig. D.2 respectively. Both correspond to the version of these models that assumes GR (namely PL5GR and PlanetGR). Full correlation plots are provided as online supplementary material.

Fig. D.1. Correlation plot of the posterior distribution function of the PL3 model using the corresponding MCMC-generated sample, and restricted to $\overline{f}$, $\overline{·}f$ as well as PL3 specific parameters.

Fig. D.2. Correlation plot of the posterior distribution function of the PlanetGR model using the corresponding MCMC-generated sample, and restricted to $\overline{f}$, $\overline{·}f$ as well as Planet specific parameters.

Appendix E: Additional material

Fig. E.1. Periodogram of the post-fit residuals, sampled at 1/T where T is the time span of observations, and ranging in frequency from 1/T to 1 day−1. Vertical dashed lines represent the best-fit planet orbital frequency, outer binary and inner binary orbital frequencies, from left to right. The vertical dotted line marks the third harmonic of the planet orbital frequency. Horizontal grey lines describe the 95% region obtained by bootstrapping white noise: the central dotted line is the median, dashed lines delineate the 2.5% and 97.5% levels, and thin dotted lines the 0.5% to 99.5% interval. For clarity, these lines have been smoothed as they are otherwise noisy. In the low-frequency part of the plot, until 5/PO ≃ 1.5 × 10−2 day−1 where PO is the outer binary period, we show the periodograms of the planet (dashed orange), Kepler (dash-dotted red), PL2/Kepl1 (dotted green), and PL3 (solid blue) models. In the higher frequency part only PL3 is shown for clarity as the difference with other models becomes much smaller.

All Tables

All Figures

Fig. 1. Comparison of numerically computed Rømer delay with first-order and Keplerian-order approximations. Approximate expressions have been least-square fitted to the numerical result. Top: Three versions as well as residuals of the least-square fit (lower panel). Bottom: Lomb-Scargle periodogram of the three versions. Vertical dashed lines mark the fundamental (black), second harmonic (grey), and third harmonic (dash-dotted grey) of the planet period. Vertical dotted lines mark the fundamental (black) and first harmonic (grey) of the outer-binary period.

In the text

	Fig. 2. Dispersion measure per time interval in the model PL3DM10. The x axis gives the time interval index, and the intervals are equal. Error bars delimit the 68% confidence region.
In the text

	Fig. 3. Galactic motion of the J0337 system during the past 500 Myr. The blue dot marks its current position. The orange dot shows the location of the Sun. The orbit was calculated with the Galactic gravitational potential and the software provided by McMillan (2017).
In the text

	Fig. 4. Velocity of the J0337 system with respect to the local galactic co-rotating frame during the last 500 Myr.
In the text

	Fig. 5. Measurements of the SEP violation parameter Δ (left column) and its absolute value |Δ| (right column) assuming the PL3 (upper row) or Planet (bottom row) models. Vertical dotted lines mark the mean value of the distributions, while the vertical dashed lines delimit the 95% confidence regions.
In the text

	Fig. 6. Timing residuals. Top: Best-fit residuals of model PL3. Bottom: Difference between residuals of best-fit PL3 and Planet models. Error bars are not shown for clarity, but the median 1-sigma uncertainty is 1.9 μs, and the mean is 2.2 μs.
In the text

Fig. 7. Reduced χ2 of the residuals between each 100-day sub-profile and the main template as a function of date. Two sub-profiles are in the time interval when imperfect polarisation calibration was performed (‘UnperfectPolCal’, see main text). The inset shows the residuals obtained for one of these sub-profiles, while others do not show discernible structure (i.e. noise). The two horizontal lines delimit the two-standard-deviation interval around the mean of the reduced χ2 distribution corresponding to σred.χ2 ≃ 0.043 (imperfect polarisation region excluded).

In the text

Fig. 8. Chaos map representing the variation of the proper mean motion |nf − ni| between two halves of a 2000-year integration on a grid of 120 × 120 initial conditions in the period-eccentricity plane of the planet. White pixels represent orbits that are so chaotic that we could not run the frequency analysis. The eccentricity grid initially went all the way to 0.9, but we show here only the 102 lower values, as nearly all pixels in the last 18 lines are white. Vertical dotted lines show the positions of resonances 8:1 to 12:1 with the outer-binary period. The cross shows the mean parameters of the planet as reported in Table E.1, and the ellipses represent the one and two-sigma confidence regions, respectively.

In the text

	Fig. 9. Evolution of the eccentricity and inclination of the planet as a function of its argument of periastron. Orbital elements are measured with respect to the orbital plane of the inner two binaries. The small dots show the trace of a 20-Myr numerical integration of the planet, as given by an N-body integrator without GR. Initial conditions are chosen near the best-fit orbit.
In the text

	Fig. 10. Probability of a planet surviving the SN explosion as a function of its pre-SN orbital period and the present-day systemic velocity of J0337. The red circle marks our default value with $v_{\mathrm{sys}}=44{\mathrm{km}\mathrm{s}}^{-1}$ and $P_{\mathrm{orb}}^{\mathrm{planet}}=1000\mathrm{days}$, yielding Pbound = 33% (see text).
In the text

	Fig. B.1. Maps representing the values of γc (left) and γs (right) as a function of (Ω1, i2) for i1 = 39deg.
In the text

	Fig. D.1. Correlation plot of the posterior distribution function of the PL3 model using the corresponding MCMC-generated sample, and restricted to $\overline{f}$, $\overline{·}f$ as well as PL3 specific parameters.
In the text

	Fig. D.2. Correlation plot of the posterior distribution function of the PlanetGR model using the corresponding MCMC-generated sample, and restricted to $\overline{f}$, $\overline{·}f$ as well as Planet specific parameters.
In the text

Fig. E.1. Periodogram of the post-fit residuals, sampled at 1/T where T is the time span of observations, and ranging in frequency from 1/T to 1 day−1. Vertical dashed lines represent the best-fit planet orbital frequency, outer binary and inner binary orbital frequencies, from left to right. The vertical dotted line marks the third harmonic of the planet orbital frequency. Horizontal grey lines describe the 95% region obtained by bootstrapping white noise: the central dotted line is the median, dashed lines delineate the 2.5% and 97.5% levels, and thin dotted lines the 0.5% to 99.5% interval. For clarity, these lines have been smoothed as they are otherwise noisy. In the low-frequency part of the plot, until 5/PO ≃ 1.5 × 10−2 day−1 where PO is the outer binary period, we show the periodograms of the planet (dashed orange), Kepler (dash-dotted red), PL2/Kepl1 (dotted green), and PL3 (solid blue) models. In the higher frequency part only PL3 is shown for clarity as the difference with other models becomes much smaller.

In the text