# Explanation of the exceptionally strong timing noise of PSR J0337+1715 by a circum-ternary planet and consequences for gravity tests

Guillaume Voisin, Ismaël Cognard, Melaine Saillenfest, Thomas Tauris, Norbert Wex, Lucas Guillemot, Gilles Theureau, P.C.C. Freire, Michael Kramer

► To cite this version:

Guillaume Voisin, Ismaël Cognard, Melaine Saillenfest, Thomas Tauris, Norbert Wex, et al. Explanation of the exceptionally strong timing noise of PSR J0337+1715 by a circum-ternary planet and consequences for gravity tests. Astronomy & Astrophysics - A&A, 2025, 693 (A143), 10.1051/0004-6361/202452100. [hal-04784066v2](https://hal.science/hal-04784066v2)

HAL Id: hal-04784066

https://hal.science/hal-04784066v2

Submitted on 23 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Explanation of the exceptionally strong timing noise of PSR J0337+1715 by a circum-ternary planet and consequences for gravity tests

G. Voisin<fnref n="1" />,*, I. Cognard<fnref n="2" />,<fnref n="3" />, M. Saillenfest<fnref n="4" />, T. M. Tauris<fnref n="5" />,<fnref n="6" />, N. Wex<fnref n="6" />, L. Guillemot<fnref n="2" />,<fnref n="3" />, G. Theureau<fnref n="2" />,<fnref n="3" />,<fnref n="1" />, P. C. C. Freire<fnref n="6" />, and M. Kramer<fnref n="6" />

Received 3 September 2024; accepted 12 November 2024

## ABSTRACT

**Context.** Timing of pulsar PSR J0337+1715 provides a unique opportunity to test the strong equivalence principle (SEP) with a strongly self-gravitating object. This is due to its unique situation in a triple stellar system with two white dwarfs.

**Aims.** Our previous study suggested the presence of a strong low-frequency residual signal in the timing data, and we set out to model this signal on a longer dataset in order to determine its nature and improve accuracy.

**Methods.** We considered three models: chromatic red noise, achromatic red noise, and a small planet in a hierarchical orbit with the triple stellar system. These models were implemented in our numerical timing model. We performed Bayesian inference of posterior distributions. Best fits were compared using information-theoretic criteria.

**Results.** We rule out chromatic red noise from dispersion-measure variations. Achromatic red noise or a planet in Keplerian orbit provide the best fits. If the residual signal is red noise, then it appears exceptionally strong. When assuming the presence of a planet, we obtained a marginal detection of mutual interactions that allowed us to constrain its mass to ~ 0.5$M_{\text{Moon}}$ as well as its inclination. The latter is intriguingly coincident with a Kozai resonance. We show that a longer observation span will ultimately lead to a clear signature of the planet model due to its mutual interactions with the triple system. We produce new limits on SEP violation: $|\Delta| < 1.5 \cdot 10^{-6}$ or $|\Delta| < 2.3 \cdot 10^{-6}$ at a 95% confidence level under the planet or red-noise hypothesis, respectively. This model dependence emphasises the need for additional data and model selection. As a by-product, we estimated a rather low supernova kick velocity of ~ 110–125km/s, strengthening the idea that it is a necessary condition for the formation of pulsar triple systems.

**Key words.** Gravitation – (Stars:) pulsars: individual PSR J0337+1715 – Stars: neutron – Radio continuum: stars – Planetary systems

## 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 extrasolar 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.2M_{\odot}$ 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.4M_{\odot}$ 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

<footnote n="1">¹ Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université de Paris, CNRS, F-92190 Meudon, France</footnote>

<footnote n="2">² Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans / CNRS, 45071 Orléans Cedex 02, France</footnote>

<footnote n="3">³ Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France</footnote>

<footnote n="4">⁴ IMCCE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Lille, 75014 Paris, France</footnote>

<footnote n="5">⁵ Dept. of Materials and Production, Aalborg University, DK-9220 Aalborg Øst, Denmark</footnote>

<footnote n="6">⁶ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote>* Email: guillaume.voisin@obspm.fr</footnote>
B1620–26 (Backer et al. 1993; Thorsett et al. 1999). The sys-
tem is composed of a likely gas giant companion orbiting a pul-
sar 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 sys-
tems: ZTF J1406+1222, which is a black-widow pulsar candi-
date that has recently been suggested to be associated with a dis-
tant low-mass stellar companion in a ~ 12, 000-year orbit (Bur-
dge et al. 2022) extending over ~ 600AU; the redback pulsar
PSR 2043+1711, which has a suspected distant stellar compan-
ion (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 gravita-
tional weak equivalence principle (GWEP; Will 2014). Together
with the local Lorentz and position invariances of gravity, it con-
stitutes 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;<fnref n="1" /> all others pre-
dict 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 vi-
olation will either rule out GR (in case of a detection) or alterna-
tive theories of gravity (in case of a non-detection).

In alternative theories of gravity, such as scalar-tensor the-
ories (Damour & Esposito-Farèse 1993; Damour & Esposito-
Farese 1996), violation of GWEP takes the form of an effective
body-dependent gravitational constant G<sub>ab</sub> = G<sub>ba</sub> = G(1 + Δ<sub>ab</sub>),
where G is the Newtonian gravitational constant and Δ<sub>ab</sub> is a vi-
olation parameter, the definition of which is theory dependent. In
weakly self-gravitating bodies this reduces to the usual discrep-
ancy 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 proportional-
ity 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<sub>ab</sub> 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 sys-
tem, 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 MI-
CROSCOPE 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 (⊕ - o - ⊙) 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. Further-
more, 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 ob-
ject 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 pul-
sar - white dwarf binary evolving in the gravitational potential of
the Galaxy. A violation of GWEP would polarise the orbital ec-
centricity of the binary towards the Galactic centre. Using one of
the most accurately timed pulsars, PSR J1713+0747, Zhu et al.
(2019) obtained |Δ<sub>pG</sub> − Δ<sub>cG</sub>| < 2 × 10<sup>−3</sup> 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 sensi-
tive 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 sys-
tem 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 in-
tegrations, Voisin et al. (2020b) independently developed a nu-
merical timing model (Nutimo) (Voisin et al. 2020b) that sup-
ports 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 nec-
essary 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 viola-
tion 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 compan-
ions. 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

<footnote n="1">¹ The only other SEP-abiding theory is Nordström’s scalar theory (Gérard 2007; Deruelle 2011); however, this was shown to be inconsistent with several experimental results on light deflection.</footnote>
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 (Niu 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 Niu 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).<fnref n="2" /> 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-1730MHz) 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 2). 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 Sec. 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

<footnote n="2">² Dataset, version of Nutimo used in this work, and results are available on Zenodo at https://doi.org/10.5281/zenodo.13899771.</footnote>
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

$$F(t_a) = \sum_{k=1}^{n} \frac{A}{k^\gamma} \sin(2k\pi vt_a + \phi_k), \quad (1)$$

where A is the amplitude of the Fourier component at the fundamental frequency $v$, $\gamma$ is the power-law index, $\phi_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 \in \{2, 3, 4, 5\}$. There are $3 + n$ fitted parameters, which are $v, A, \gamma, \phi_{k \in \{1, n\}}$. We remark that, contrary to the stochastic process approach, we fit for the fundamental frequency $v$. 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 $\approx 299$ days, which is much smaller than the residual signal timescale of $\sim 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 \approx 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_\pi$ orbiting the triple system in a hierarchical orbit, that is with a period $P_\pi \gg 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_\pi \ll M_\odot$. The only measurable effect given the amplitude of the signal (a few $\mu$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 $(\mathbf{R}, \mathbf{V})_{k \in [p, i, o, \pi]}$ be the positions and velocities of the pulsar (p), inner white dwarf (i), outer white dwarf (o), and planet ($\pi$) 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 $\Pi = \{t, \pi\}$, the centre of mass of which is also centre of mass of the whole system.

One can decompose the position of the pulsar as

$$ \mathbf{R}_p = \mathbf{R}_{p/b} + \mathbf{R}_{b/t} + \mathbf{R}_t, \quad (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{pmatrix} (\mathbf{R}, \mathbf{V})_{p/b} \\ (\mathbf{R}, \mathbf{V})_{b/t} \\ (\mathbf{R}, \mathbf{V})_t \\ m_p, m_i, m_o, m_\pi \end{pmatrix} \leftrightarrow \begin{pmatrix} a_p, e_l, \omega_l, t_{ascp}, i_l, \Omega_b \\ a_b, e_O, \omega_O, t_{ascb}, i_O, \Omega_b \\ a_t, e_{\Pi}, \omega_{\Pi}, t_{asct}, i_{\Pi}, \Omega_t \\ m_i/m_p, P_I, P_O, P_{\Pi} \end{pmatrix}, \quad (3) $$

where $a_x, e_x, \omega_x, t_{ascx}, i_x, \Omega_x$ are the semi-major axis, orbital eccentricity, argument of periastron, time of passage at the ascending node,<fnref n="3" /> 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 \in \{\text{p, b, t}\}$ and binaries $X \in \{\text{I, O, \Pi}\}$. 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,\pi}$ are related to the mass ratio $m_i/m_p$ and orbital periods $P_{I,O,\Pi}$ 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_X \cos \omega_X, e_X \sin \omega_X$ particularly adapted to low-eccentricities (e.g. Lange et al. 2001) and the use of $\delta i = i_l - i_O, \delta \Omega = \Omega_b - \Omega_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,

$$ \Delta_R = -\frac{n_{\odot} \cdot \mathbf{R}_p}{c} = -\frac{1}{c} (\mathbf{n}_{\odot} \cdot \mathbf{R}_{p/b} + \mathbf{n}_{\odot} \cdot \mathbf{R}_{b/t} + \mathbf{n}_{\odot} \cdot \mathbf{R}_t), \quad (4) $$

where $\mathbf{n}_{\odot}$ 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, $\mathbf{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_X \sin i_X$, the orbital eccentricity $e_X$ and argument of periastron $\omega_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:

$$ a_i^3 n_{\Pi}^2 = G \frac{m_{\pi}^3}{(m_t + m_{\pi})^2}, \quad (5) $$

where $n_{\Pi} = 2\pi/P_{\Pi}$ 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 $\sin i_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 $\mathbf{R}_t = \mathbf{R}_t^K + \delta \mathbf{R}_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 $\delta \Delta_R = -\mathbf{n}_{\odot} \cdot \delta \mathbf{R}_t/c$ of the form, Eq. (B.22),

$$ \delta \Delta_R = \alpha t (\gamma_c \cos n_O t + \gamma_s \sin n_O t), \quad (6) $$

where $n_O = 2\pi/P_O$ is the mean motion of the outer binary and $\alpha, \gamma_c, \gamma_s$ are constants defined in (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

<footnote n="3">³ For practical purposes we use only an approximation of the time of passage at the ascending node, which is defined relative to the time of periastron passage $t_p$ by $t_{ascx} = t_{px} - P_X \omega_x / 2\pi$.</footnote>
<table><caption>**Table 1.** Mean values of model-specific parameters of the MCMC runs assuming the PL3 (left) and Planet (right) models with their 68% confidence intervals.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>PL3GR</th><th>PlanetGR</th><th>Symbol</th><th>Parameter</th></tr></thead><tbody><tr><td colspan="6" style="text-align:center; font-weight:bold;">Fitted values</td></tr><tr><td>Fundamental frequency (day<sup>-1</sup>)</td><td>v</td><td>3.44(28)<sup>+61</sup><sub>-54</sub> × 10<sup>-4</sup></td><td>3310<sup>+280</sup><sub>-260</sub></td><td>P<sub>II</sub></td><td>Orbital period (d)</td></tr><tr><td>Amplitude at v (µs)</td><td>A<sub>v</sub></td><td>4.4(82)<sup>+22</sup><sub>-19</sub></td><td>(6.5)<sup>+1.6</sup><sub>-1.4</sub> × 10<sup>-6</sup></td><td>a<sub>t</sub> sin i<sub>II</sub></td><td>Projected semi-major axis (lt-s)</td></tr><tr><td>Power-law index</td><td>γ</td><td>2.73(26)<sup>+68</sup><sub>-61</sub></td><td>-(1.5)<sup>+1.4</sup><sub>-1.6</sub> × 10<sup>-5</sup></td><td>a<sub>t</sub> cos i<sub>II</sub></td><td>Co-projected semi-major axis (lt-s)</td></tr><tr><td>Phase 1 (rad)</td><td>φ<sub>1</sub></td><td>-0.21(54)<sup>+19</sup><sub>-17</sub></td><td>0.2(51)<sup>+32</sup><sub>-37</sub></td><td>e<sub>II</sub> sin ω<sub>t</sub></td><td>Laplace-Lagrange</td></tr><tr><td>Phase 2 (rad)</td><td>φ<sub>2</sub></td><td>-0.77(74)<sup>+68</sup><sub>-61</sub></td><td>(5.7)<sup>+2.4</sup><sub>-2.4</sub> × 10<sup>-2</sup></td><td>e<sub>II</sub> cos ω<sub>t</sub></td><td>Laplace-Lagrange</td></tr><tr><td>Phase 3 (rad)</td><td>φ<sub>3</sub></td><td>-1.42(21)<sup>+68</sup><sub>-61</sub></td><td>56549<sup>+3333</sup><sub>-20</sub></td><td>t<sub>asct</sub></td><td>Time of ascending node (MJD)</td></tr><tr><td></td><td></td><td></td><td>1.(25)<sup>+62</sup><sub>-53</sub> × 10<sup>2</sup></td><td>Ω<sub>t</sub></td><td>Longitude of planet ascending node (°)</td></tr><tr><td colspan="6" style="text-align:center; font-weight:bold;">Derived values</td></tr><tr><td>Fundamental period (days)</td><td>1/v</td><td>2904.6<sup>+4.6</sup><sub>-5.1</sub></td><td>(1.65)<sup>+1.6</sup><sub>-0.89</sub> × 10<sup>-5</sup></td><td>a<sub>π</sub></td><td>Semi-major axis (lt-s)</td></tr><tr><td>Amplitude at 1yr<sup>-1</sup> (µs)</td><td>A<sub>1yr<sup>-1</sup></sub></td><td>1.55(16)<sup>+61</sup><sub>-68</sub> × 10<sup>-2</sup></td><td>157<sup>+11</sup><sub>-54</sub></td><td>i<sub>II</sub></td><td>Orbital inclination w.r.t. plane of sky (°)</td></tr><tr><td></td><td></td><td></td><td>119<sup>+16</sup><sub>-42</sub> × 10<sup>2</sup></td><td>δi<sub>II</sub></td><td>Orbital inclination w.r.t. triple system (°)</td></tr><tr><td></td><td></td><td></td><td>0.2(57)<sup>+29</sup><sub>-31</sub></td><td>e<sub>II</sub></td><td>Orbital eccentricity</td></tr><tr><td></td><td></td><td></td><td>77.2<sup>+6.0</sup><sub>-7.2</sub></td><td>ω<sub>t</sub></td><td>Longitude of periastron (°)</td></tr><tr><td></td><td></td><td></td><td>57258<sup>+3270</sup><sub>-77</sub> × 10<sup>4</sup></td><td>t<sub>pII</sub></td><td>Time of periastron passage (MJD)</td></tr><tr><td></td><td></td><td></td><td>(1.23)<sup>+1.1</sup><sub>-0.66</sub> × 10<sup>-8</sup></td><td>m<sub>π</sub></td><td>Planet-companion mass (M<sub>☉</sub>)</td></tr></tbody></table>

**Notes.** Uncertainties apply to the digits between parenthesis and delimit the 68% median confidence region, namely the interval excluding 16% of the distribution above and below it. The central value quoted is the best-fit value. The upper-case index II refers to the planet binary. The lower-case index π refers to the planet, and t refers to the triple-system centre of mass. Data spans MJD 56492 – 59480. Solar System ephemeris is DE430. Parameters common to both models are reported in Table E.1.

degenerate with other components of the timing formula and can in principle allow for the separate measurement of $\alpha\gamma_c$ and $\alpha\gamma_s$.

Importantly, $\gamma_c$ and $\gamma_s$ are independent functions of inclination $i_{II}$ and longitude of ascending node $\Omega_t$, while $\alpha$ 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 + \Delta_{ab})$, where G is the Newtonian gravitational constant (see also Sect. 1 and Eq. (20) of Voisin et al. (2020b)). In GR, $\Delta_{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 $\Delta \equiv \Delta_{pb}$, where $b \in \{i, o, \pi\}$.

## 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; Kapan 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, $\hat{\sigma}_{\dot{m}}$ and $\hat{\sigma}_{\dot{\delta}}$ 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 $\sigma$, $\hat{\sigma}_{\dot{m}}/\sigma$ and $\hat{\sigma}_{\dot{\delta}}/\sigma$ respectively, are used as 'convergence ratios' (Dunkley et al. 2005)$^4$. 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 Criterion (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

<footnote>$^4$ We note that (Dunkley et al. 2005) estimate directly the variance of the mean of the whole chain using a more refined spectral method while we estimate variances of a sample of sub-chains by the usual variance estimator. These sub-chain variance estimates ought to be larger than full-chain variances since the estimator scales as $\sim 1/n$ where n is sample size. Therefore, albeit simpler, this method is conservative.</footnote>
<table><caption>**Table 2.** Best-fit statistical properties of the models fitted to the dataset of this paper.</caption><thead><tr><th>Model</th><th>Npar</th><th>χ<sup>2</sup></th><th>R<sub>χ</sub><sup>2</sup></th><th>ΔAIC</th><th>ΔBIC</th></tr></thead><tbody><tr><td>Planet</td><td>32</td><td>+2.4</td><td>+0.0003</td><td>+4.4</td><td>+12</td></tr><tr><td>Kepler</td><td>30</td><td>+4.0</td><td>0.0002</td><td>+2.0</td><td>-5.4</td></tr><tr><td>PL2/Kep11</td><td>30</td><td>+19</td><td>+0.0014</td><td>+17</td><td>+9.6</td></tr><tr><td>PL3</td><td>31</td><td>15511.9</td><td>1.2466</td><td>15573.9</td><td>15804.3</td></tr><tr><td>PL4</td><td>32</td><td>+4.8</td><td>+0.0005</td><td>+6.8</td><td>+14</td></tr><tr><td>PL5</td><td>33</td><td>-5.4</td><td>-0.0002</td><td>-1.4</td><td>+14</td></tr><tr><td>PL3DM10</td><td>39</td><td>-17</td><td>-0.0006</td><td>-1.3</td><td>+58</td></tr></tbody></table>

**Notes.** Npar is the number of parameters of the model, $\chi^2$ the minimum chi squared, $R_{\chi}^2$ is equal to $\chi^2/N_{\text{dof}}$ where $N_{\text{dof}} = 12474 - \text{Npar}$. AIC and BIC are the Akaike and Bayesian information criteria, respectively. The ‘PL3’ line is the reference fit subtracted from all the other lines. All models assume GR. Models are described in Sec. 4.1: ‘Planet’ is the full 4-body numerical model; ‘Kepler’ is the purely Keplerian model without mutual interactions; ‘PLn’ are models with n sinusoidal harmonics following Eq. (1) in order to account for red noise; PL2 is equivalent to a planet in Keplerian orbit with the centre of mass of the triple system at leading order in eccentricity (no mutual interactions) and thus denoted ‘PL2/Kep11’; PL3DM10 fits separate DM within ten evenly distributed intervals on top of three harmonics (We note the other models have DM and DM’ parameters, which have been removed here).

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’)<fnref n="5" /> 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/Kep11’. Given the small difference between Kepler and PL2/Kep11, no additional MCMC was done for Kepler, and its best-fitting solutions are based on a Minuit refit of PL2/Kep11, thus saving on computing power.

Assuming the presence of a planet, the PL2/Kep11 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 \sim 0.25$ and a projected semi-major axis $x \sim 4\mu\text{s}$ (see Table 1), the amplitude of this harmonic is only about $\circ(xe^2) \sim 0.25\mu\text{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/Kep11 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/Kep11 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.

The Planet model has a moderately better $\chi^2$ 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/Kep11 concerning BIC. This suggests that mutual interactions may be detected, since the $\chi^2$ is better, but only marginally, which is in line with the large uncertainty on the planet mass (see Table 1).

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

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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 ($\Delta\text{AIC} \approx -1.38$). Thus, we make it our reference achromatic red-noise model in the following. PL3 has the best $\chi^2$ of all models, moderately better than Planet ($\Delta\chi^2 \approx -2.4$), and a moderately better AIC than Kepler ($\Delta\text{AIC} \approx -2.0$) but significantly worse BIC ($\Delta\text{BIC} \approx +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 $\chi^2$ 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.

We also ran a model called PL3DM10, for which the MCMC only sampled the six achromatic red-noise parameters, the ten DMx DM parameters as well as the pulsar spin parameters $\tilde{f}, \tilde{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 DMx. 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 DMx values is compa-

<footnote n="5">⁵ We use ‘Planet’ with capital P whenever referring to the model.</footnote>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

rable 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 2.

The timing signal introduced by both models are shown in Fig. 1. In that figure we have also represented the signal from the PL2/Kep11 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/Kep11, 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 $\delta i$ and relative longitude of ascending nodes $\delta\Omega$ being consistent with zero within error bars in both models. We also note that the tension in declination proper motion $\mu_{\delta}$ compared to its GAIA prior that was observed in Voisin et al. (2020b) has disappeared in the present results.

The longitude of ascending node $\Omega_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 $\Omega_b$ by timing.<fnref n="6" /> However, the MCMC seemed unable to find a solution accounting for this peak while staying in the vicinity of the initial value of $\Omega_b$, and the assumption was made of a local minimum around the $\Omega_b$ value from Voisin et al. (2020b). The MCMC chain was restarted with values of $\Omega_b$ offset by ±90 deg and 180 deg 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 180 deg 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 ±180 deg 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

<footnote n="6">⁶ It is also possible to use scintillation arcs (Reardon et al. 2020).</footnote>
in the range ~ [30, 55] km s⁻¹ in the entire integration span (see Fig. 4). This finding is of particular importance for Sec. 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.

<figure><img src="image_5.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Velocity of the J0337 system with respect to the local galactic co-rotating frame during the last 500 Myr.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

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

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

$$ \Delta = 0.97^{+1.58}_{-1.48} \times 10^{-6} \quad \text{and} \quad |\Delta| < 2.29 \times 10^{-6} \text{ (PL3)}, \quad (7) $$

$$ \Delta = 0.51^{+1.14}_{-1.28} \times 10^{-6} \quad \text{and} \quad |\Delta| < 1.46 \times 10^{-6} \text{ (Planet).} \quad (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,
<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

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.

## 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<fnref n="7" /> 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}^2(ν/γ_{GP})^{-γ_{GP}}$ where $α = 1/12π^2$ 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 $ν = 1yr^{-1}$ 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} \sim A(T/α)^{1/2} (ν/γ_{GP})^{γ_{GP}/2}$ where we used the notations of Eq. (1), T is the time span of observations,

and $γ_{GP} \sim 2γ$. Using the results of our PL3 model, we obtain<fnref n="8" /> $A_{GP} \sim α^{-1/2}0.04µs\,yr^{1/2}$ and $γ_{GP} \sim 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@8yr), is only surpassed by PSR J1824-2452A with ~8µs@8yr and $γ_{GP} \approx 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@8yr and $γ_{GP} \approx 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} erg/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}Hz/s$ and $\dot{f} \sim -4 \times 10^{-14}Hz/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 $α \approx -1.4$; $β \approx 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 $τ = |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} erg/s$, spin-down rate of $\dot{f} \sim -2 \times 10^{-15}Hz/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 $χ^2$, as shown in Fig. 7. The re-

<footnote n="7">⁷ The six EPTA pulsars in the cited reference are also among the NANOGrav 47 pulsars, and we count 6 out of 26 pulsars in common with the PPTA dataset.</footnote>

<footnote n="8">⁸ Equivalently $\log_{10}(A_{GP}/yr^{3/2}) \sim -13.9$ with $α = (12π^2)^{-1}$ or $\log_{10}(A_{GP}/yr^{3/2}) \sim -14.9$ with $α = 1$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Reduced $\chi^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 $\chi^2$ distribution corresponding to $\sigma_{\text{red},\chi^2} \approx 0.043$ (imperfect polarisation region excluded).</figcaption></figure>

duced $\chi^2$ 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 $\chi^2$, we show their two-standard-deviation interval in Fig. 7 computed on the clean $\chi^2$ sample only. All but one element fall into that region, in line with expectations. Indeed, each profile possesses ~2000 degrees of freedom, thus the $\chi^2$ distribution is well approximated by a Gaussian law and the above-mentioned interval corresponds to a probability of 95%.

## 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 $\chi^2$ and AIC. The full Planet model does not improve sufficiently its $\chi^2$ 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 \times 10^{-5} M_{\text{Jup.}}$, ~0.004$M_{\text{Earth}}$ or ~0.4$M_{\text{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<fnref n="9" />.

We have evaluated the short-term stability of the orbit of the candidate planet by applying the frequency analysis method of Laskar (1988, 1990, 1993, 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.

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^{-5}$). 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_{\Pi} = 119_{-42}^{+16}^{\circ}$ (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,<fnref n="10" /> whose location can be approximated by $\cos(\delta i_{\Pi}) = -1/\sqrt{5}$, that is, $\delta i_{\Pi} \approx 116.6^{\circ}$ (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 $\delta i_{\Pi}$. As

<footnote n="9">⁹ The catalogue at exoplanet.eu reports one lighter planet around the white dwarf WD 1145+017. However, it seems that this object could be composed of multiple disintegrating planetesimals Croll et al. (2017).</footnote>

<footnote n="10">¹⁰ This is the usually called Lidov-Kozai effect. We follow here Ito & Ohtsuka (2024) who showed that earlier work by von Zeipel should also be acknowledged. See references in the text.</footnote>
<figure><img src="image_10.png" /></figure>

**Fig. 8.** Chaos map representing the variation of the proper mean motion $|n_f - n_i|$ 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.

the planet is in the Kozai resonance, its eccentricity and inclination are affected by correlated oscillations, and its argument of periastron oscillates around $\pi/2$ (instead of circulation between 0 and $2\pi$ 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.

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 $\omega_\pi$ it implies. Indeed, an argument of periastron confined near $\pi/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

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

**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.

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_{\text{planet}} = 1.23 \times 10^{-8} M_{\odot}$. At stage 4, just prior to the SN explosion, we can calculate the probability that this planet survived the kinematic effects of the SN.

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):
<figure><img src="image_12.png" /><figcaption>**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_{sys} = 44$ km s⁻¹ and $P_{orb}^{planet} = 1000$ days, yielding $P_{bound} = 33\%$ (see text).</figcaption></figure>

$$P_{\text{bound}} = \frac{1}{2} \left\{ 1 + \left[ \frac{1 - 2\Delta M/M - (v_{\text{sys}}/v_{\text{rel}})^2}{2(v_{\text{sys}}/v_{\text{rel}})} \right] \right\}, \quad (9)$$

where M is the total mass of the system, ΔM is the amount of material ejected in the SN, $v_{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_{\odot}$ (1.70 + 1.10 + 1.30 $M_{\odot}$), $\Delta M = 0.42 M_{\odot}$, $v_{sys} = 44$ km s⁻¹ and $P_{orb}^{planet} = 1000$ days. The latter two values are our default values: $v_{sys} = 44$ km s⁻¹ 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_{orb}^{planet} = 1000$ days is selected as a typical value from our range of solutions to the SN event based on Monte Carlo simulations with $2 \times 10^6$ 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_0 = M_0/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_{\Pi} = 119_{-42}^{+160}^{\circ}$ 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ᵢ+WD₀ and “body 2”: planet; step 2 consisting of “body 1”: NS+WDᵢ and “body 2”: WD₀; step 3 consisting of “body 1”: NS and “body 2”: WDᵢ) we find a rough constraint on the NS kick of $w \approx 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 \approx 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 2 for software and data release.)

We have shown in Sec. 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 Sec.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 Tabs. 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; Nițu 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 Sec. 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 Sec. 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–125km/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 Sec. 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 Sec. 4.4, we have shown that the limit is somewhat model dependent, with $|Δ| < 2.29 × 10^{-6}$ assuming the best red-noise model and $|Δ| < 1.46 × 10^{-6}$ 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.

**Acknowledgements.** 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

Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716

Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021, ApJ Supplement Series, 252, 4, aDS Bibcode: 2021ApJS.252....4A

Archibald, A. M., Gusinskaia, N. V., Hessels, J. W. T., et al. 2018, Nature, 559, 73

Backer, D. C., Foster, R. S., & Sallmen, S. 1993, Nature, 365, 817, aDS Bibcode: 1993Natur.365..817B

Bailes, M., Bates, S. D., Bhalerao, V., et al. 2011, Science, 333, 1717, aDS Bibcode: 2011Sci...333.1717B

Bassa, C. G., Janssen, G. H., Stappers, B. W., et al. 2016, MNRAS, 460, 2207, publisher: OUP ADS Bibcode: 2016MNRAS.460.2207B

Bear, E. & Soker, N. 2014, MNRAS, 444, 1698

Behrens, E. A., Ransom, S. M., Madison, D. R., et al. 2020, ApJ Letters, 893, L8

Beuermann, K., Hessman, F. V., Dreizler, S., et al. 2010, A&A, 521, L60

Beutler, G. 2004, Methods of Celestial Mechanics: Volume I: Physical, Mathematical, and Numerical Principles, 1st edn. (Berlin; New York: Springer)

Blandford, R. & Teukolsky, S. A. 1976, ApJ, 205, 580

Boyles, J., Lynch, R. S., Ransom, S. M., et al. 2013, ApJ, 763, 80, aDS Bibcode: 2013ApJ...763....80B

Burdge, K. B., Marsh, T. R., Fuller, J., et al. 2022, Nature, 605, 41, aDS Bibcode: 2022Natur.605...41B

Chalumeau, A., Babak, S., Petiteau, A., et al. 2022, MNRAS, 509, 5538, aDS Bibcode: 2022MNRAS.509.5538C

Champion, D. J., Ransom, S. M., Lazarus, P., et al. 2008, Science, 320, 1309

Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970
Cognard, I., Theureau, G., Guillemot, L., et al. 2013, in SF2A-2013: Proceedings of the Annual meeting of the French Society of A&A, ed. L. Cambresy, F. Martins, E. Nuss, & A. Palacios, 327–330

Croll, B., Dalba, P. A., Vanderburg, A., et al. 2017, ApJ, 836, 82, aDS Bibcode: 2017ApJ...836...82C

Damour, T. & Deruelle, N. 1985, Ann. Inst. Henri Poincaré Phys. Théor., Vol. 43, No. 1, p. 107 - 132, 43, 107

Damour, T. & Esposito-Farèse, G. 1993, Phys. Rev. Lett., 70, 2220, publisher: American Physical Society

Damour, T. & Esposito-Farese, G. 1996, Phys. Rev. D, 54, 1474, arXiv: gr-qc/9602056

Damour, T. & Schäfer, G. 1991, Phys. Rev. Lett., 66, 2549, publisher: American Physical Society

Damour, T. & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840

Deruelle, N. 2011, Gen. Relativ. Gravitation, 43, 3337

Desvignes, G., Barott, W. C., Cognard, I., Lespagnol, P., & Theureau, G. 2011, in American Institute of Physics Conference Series, Vol. 1357, Radio Pulsars: An Astrophysical Key to Unlock the Secrets of the Universe, ed. M. Burgay, N. D'Amico, P. Esposito, A. Pellizzoni, & A. Possenti, 349–350

Developers, T. S. 2022, sagemath/sage: 9.5

Donlon, II, T., Chakrabarti, S., Lam, M. T., et al. 2024, The Anomalous Acceleration of PSR J2043+1711: Long-Period Orbital Companion or Stellar Flyby?, publication Title: arXiv e-prints ADS Bibcode: 2024arXiv240706482D

Dražkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in Astron Soc Pac Conf Ser, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 717

Dunkley, J., Bucher, M., Ferreira, P. G., Moodley, K., & Skordis, C. 2005, MNRAS, 356, 925

Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJ, 915, L12

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publ. Astron. Soc. Pac., 125, 306

Freire, P. C. C., Bassa, C. G., Wex, N., et al. 2011, MNRAS, 412, 2763

Freire, P. C. C., Kramer, M., & Wex, N. 2012, Classical Quantum Gravity, 29, 184007

Freire, P. C. C. & Wex, N. 2024, Living Rev. Relativ., 27, 5

Gallardo, T., Hugo, G., & Pais, P. 2012, Icarus, 220, 392, aDS Bibcode: 2012Icar.220...392G

Gastineau, M. & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194

Gérard, J. M. 2007, Classical Quantum Gravity, 24, 1867

Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2021, MNRAS, 502, 478, aDS Bibcode: 2021MNRAS.502...478G

Goodman, J. & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Grunthal, K., Krishnan, V. V., Freire, P. C. C., et al. 2024, Triple trouble with PSR J1618-3921: Mass measurements and orbital dynamics of an eccentric millisecond pulsar, arXiv:2409.03615 [astro-ph]

Guillemot, L., Cognard, I., van Straten, W., Theureau, G., & Gérard, J. A&A, 678, A79, arXiv:2308.01599 [astro-ph]

Hills, J. G. 1983, ApJ, 267, 322

Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974

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

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

Hofmann, F. & Müller, J. 2018, Classical Quantum Gravity, 35, 035015, publisher: IOP Publishing

Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, Publications of the Astronomical Society of Australia, 21, 302, aDS Bibcode: 2004PASA...21...302H

Ito, T. & Ohtsuka, K. 2024, The Lidov-Kozai Oscillation and Hugo von Zeipel, pages: 21193492 Bytes arXiv:1911.03984 [astro-ph]

James, F. & Roos, M. 1975, Computer Physics Communications, 10, 343

Kaplan, D. L., Kupfer, T., Nice, D. J., et al. 2016, ApJ, 826, 86, publisher: IOP ADS Bibcode: 2016ApJ...826...86K

Kaplan, D. L., van Kerkwijk, M. H., Koester, D., et al. 2014, ApJ Letters, 783, L23

Kerr, M., Johnston, S., Hobbs, G., & Shannon, R. M. 2015, ApJ, 809, L11, aDS Bibcode: 2015ApJ...809L...11K

Knight, H. S., Bailes, M., Manchester, R. N., & Ord, S. M. 2006, ApJ, 653, 580, aDS Bibcode: 2006ApJ...653...580K

Konacki, M. & Wolszczan, A. 2003, ApJ Letters, 591, L147

Kopeikin, S. M. 1995, ApJ Letters, 439, L5

Kozai, Y. 1962, The Astronomical Journal, 67, 591, aDS Bibcode: 1962AJ...67...591K

Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2021, Phys. Rev. X, 11, 041050

Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274

Laskar, J. 1988, A&A, 198, 341

Laskar, J. 1990, Icarus, 88, 266, aDS Bibcode: 1990Icar...88...266L

Laskar, J. 1993, Physica D, 67, 257

Laskar, J. 2005, in Hamiltonian Systems and Fourier Analysis: New Prospects for Gravitational Dynamics, ed. D. Benest, C. Froeschle, & E. Lega (Cambridge Scientific Publishers)

Lidov, M. L. 1962, Planetary and Space Science, 9, 719

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2, publisher: EDP Sciences

Liu, Z.-W., Tauris, T. M., Röpke, F. K., et al. 2015, A&A, 584, A11

Lynch, R. S., Boyles, J., Ransom, S. M., et al. 2013, ApJ, 763, 81, aDS Bibcode: 2013ApJ...763...81L

Lyne, A. & Graham-Smith, F. 2012, Pulsar Astronomy, 4th edn. (Cambridge : New York: Cambridge University Press)

Lyne, A., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science, 329, 408

Manchester, R. N. 2017, J. Astrophys. Astron., 38, 42

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, The Astronomical Journal, 129, 1993

McMillan, P. J. 2017, MNRAS, 465, 76

Melatos, A. & Link, B. 2014, MNRAS, 437, 21

Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics, publication Title: Solar system dynamics by C.D. Murray and S.F. McDermott. (Cambridge ADS Bibcode: 1999ssd.book......M)

Nieder, L., Kerr, M., Clark, C. J., et al. 2022, ApJ, 931, L3, aDS Bibcode: 2022ApJ...931L...3N

Nitu, I. C., Keith, M. J., Stappers, B. W., Lyne, A. G., & Mickaliger, M. B. 2022, MNRAS, 512, 2446

Nordtvedt, K. 1968, Physical Review, 170, 1186

Pijloo, J. T., Caputo, D. P., & Portegies Zwart, S. F. 2012, MNRAS, 424, 2914

Portegies Zwart, S., van den Heuvel, E. P. J., van Leeuwen, J., & Nelemans, G. 2011, ApJ, 734, 55

Ransom, S. M., Stairs, I. H., Archibald, A. M., et al. 2014, Nature, 505, 520

Reardon, D. J., Coles, W. A., Bailes, M., et al. 2020, ApJ, 904, 104, aDS Bibcode: 2020ApJ...904...104R

Reardon, D. J., Shannon, R. M., Cameron, A. D., et al. 2021, MNRAS, 507, 2137, aDS Bibcode: 2021MNRAS.507.2137R

Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27, aDS Bibcode: 2021ApJ...918L...27R

Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2016, Celest. Mech. Dyn. Astron., 126, 369, aDS Bibcode: 2016CeMDA.126...369S

Schleicher, D. R. G. & Dreizler, S. 2014, A&A, 563, A61

Schwarz, G. 1978, Ann. Stat., 6, 461, publisher: Institute of Mathematical Statistics

Shannon, R. M. & Cordes, J. M. 2010, ApJ, 725, 1607

Shannon, R. M., Cordes, J. M., Metcalfe, T. S., et al. 2013, ApJ, 766, 5

Shao, L. 2016, Phys. Rev. D, 93, 084023, arXiv: 1602.05725

Sigurdsson, S., Richer, H. B., Hansen, B. M., Stairs, I. H., & Thorsett, S. E. 2003, Science, 301, 193, publisher: American Association for the Advancement of Science

Sigurdsson, S. & Thorsett, S. E. 2005, in Astron Soc Pac Conf Ser, Vol. 328, Binary Radio Pulsars, ed. F. A. Rasio & I. H. Stairs, 213

Susobhanan, A., Gopakumar, A., Joshi, B. C., & Kumar, R. 2018, MNRAS, 480, 5260

Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170

Tauris, T. M. & van den Heuvel, E. P. J. 2014, ApJ Letters, 781, L13

Tauris, T. M. & van den Heuvel, E. P. J. 2023, Physics of Binary Star Evolution. From Stars to X-ray Binaries and Gravitational Wave Sources (Princeton University Press)

Thorsett, S. E., Arzoumanian, Z., Camilo, F., & Lyne, A. G. 1999, ApJ, 523, 763, publisher: IOP ADS Bibcode: 1999ApJ...523...763T

Touboul, P., Métris, G., Rodrigues, M., et al. 2019, Classical Quantum Gravity, 36, 225006

Touboul, P., Métris, G., Rodrigues, M., et al. 2022, Phys. Rev. Lett., 129, 121102

Tsang, D. & Gourgouliatos, K. N. 2013, ApJ, 773, L17, publisher: American Astronomical Society

van Straten, W. 2006, ApJ, 642, 1004

Vleeschower, L., Corongiu, A., Stappers, B. W., et al. 2024, MNRAS, 530, 1436

Voisin, G. 2017, Theses, Université de recherche Paris Sciences et Lettres, https://hal.archives-ouvertes.fr/tel-01677325

Voisin, G., Breton, R. P., & Summers, C. 2020a, MNRAS, 492, 1550

Voisin, G., Cognard, I., Freire, P., et al. 2022, One pulsar, two white dwarfs, and a planet confirming the strong equivalence principle, arXiv:2205.09345

Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020b, A&A, 638, A24

von Zeipel, H. 1909, Astronomische Nachrichten, 183, 345, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/asna.19091832202

Wheeler, J. C., Lecar, M., & McKee, C. F. 1975, ApJ, 200, 145

Will, C. M. 2014, Living Rev. Relativ., 17, arXiv: 1403.7377

Wolszczan, A. & Frail, D. A. 1992, Nature, 355, 145

Zhu, W. W., Desvignes, G., Wex, N., et al. 2019, MNRAS, 482, 3249

Zhu, W. W., Stairs, I. H., Demorest, P. B., et al. 2015, ApJ, 809, 41, aDS Bibcode: 2015ApJ...809...41Z
Appendix A: Second-order Keplerian Røemer delay

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

$$
\Delta_R = x(\sin \omega (\cos E - e) + \sqrt{1 - e^2} \cos \omega \sin E), \quad (A.1)
$$

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

$$
E - e \sin E = \sigma, \tag{A.2}
$$

where $\sigma = n(t - T_p)$ with time $t$, time of passage at periastron $T_p$
and $n = 2\pi/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} = \sigma + e \sin E_n$. In-
jecting $E_2$ into Eq. A.1 and Taylor expanding, we obtained the
second-order Røemer delay,

$$
\begin{equation}
\begin{split}
\Delta_R = {}& x[\sin \phi \\
             & \quad + \frac{e}{2}(\cos \omega \sin 2\phi - \sin \omega \cos 2\phi - 3 \sin \omega) \\
             & \quad - \frac{e^2}{8} \{ (\cos 2\omega + 4) \sin \phi - \sin 2\omega \cos \phi \\
             & \quad \phantom{\quad - \frac{e^2}{8} \{ (} - 3 \cos 2\omega \sin 3\phi + 3 \sin 2\omega \cos 3\phi \} ],
\end{split}
\tag{A.3}
\end{equation}
$$

where we define $\phi = n(t - T_{\text{asc}})$ with $T_{\text{asc}} = T_p - \omega/n$, which corresponds to the time of ascending node at order $e^0$. The first two lines correspond to the small-eccentricity timing model of Lange et al. (2001) where one can identify the Laplace-Lagrange parameters $e \sin \omega$ and $e \cos \omega$. In all rigour we left here the term in $\sim e \cos \omega$, 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 ac-
count second-order terms is i) to modify the coefficients in front
of first-harmonic terms and ii) to add signal at the third har-
monic frequency. Without the third-harmonic contribution, i)
would merely amount to an effective eccentricity and argument
of periastron. However thanks to third harmonics this degen-
eracy 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<sub>p</sub>, n
can be measured independently at this order up to an error
⊙(e<sup>2</sup>). The equivalence between PL2 as given by Eq. (1), and
this model approximately goes as follows:

$$
n = 2\pi\nu, \tag{A.4}
$$

$$
x = A, \tag{A.5}
$$

$$
T_p = T_{\text{ref}} - \frac{\phi_2 - \phi_1}{\nu}, \qquad (\text{A.6})
$$

$$
\omega = 2\phi_1 - \phi_2, \tag{A.7}
$$

$$
e = 2^{1-\gamma}. \tag{A.8}
$$

The approximation lies in the fact that we neglected the differ-
ence between emission time and arrival time, Δ = t_a − t_e Einstein
delay excluded, such that Eqs. (A.4)-(A.8) are valid up to cor-
rections 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 re-
quires 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 corre-
sponds to the inner binary, the second one to the outer white
dwarf and the last one to the planet. The extent of the inner bi-
nary 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 $\mathbf{R}_{0,1,2}$, $\mathbf{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 ($\mathbf{R}_0$, $\mathbf{V}_0$) 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),

$$
r_0 = \mathbf{B}_2, \tag{B.1}
$$

$$
r_i = R_i - B_{i-1} \quad (i > 0), \tag{B.2}
$$

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

$$
\mathbf{B}_i = \frac{1}{\sigma_i} \sum_{k=0}^{i} M_k \mathbf{R}_k, \qquad (\text{B.3})
$$

$$
\text{with } \sigma_i = \sum_{k=0}^{i} M_k.
$$

With these coordinates $\mathbf{r}_1$ is the vector connecting the inner
binary to the outer white dwarf, which we call the outer binary,
and $\mathbf{r}_2$ 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:

$$
\epsilon \sim M_2 / M_{0,1} \ll 1, \qquad (B.4)
$$

$$
\eta \sim r_2 / r_1 \sim R_2 / R_{0,1} \ll 1, \quad (B.5)
$$

where we use the notation $V = \|V\|$.

Using canonical coordinates and their conjugate momenta
**P**<sub>0,1,2</sub>, the Hamiltonian of the system is given by

$$
H = \frac{1}{2} \sum_{k=0}^{2} \frac{\mathbf{P}_k^2}{M_k} - \sum_{0 \le k < l \le 2} \frac{G M_k M_l}{\|\mathbf{R}_l - \mathbf{R}_k\|}, \quad (B.6)
$$

where G is the gravitational constant.

We now decompose the Hamiltonian in an dominant inte-
grable part $H_0$ and a perturbing term $H_1$ 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 sub-
stitution $\mathbf{P}_k \rightarrow \mathbf{p}_k$, $M_k \rightarrow m_k$ where $m_0 = \sigma_2$ (the total mass
of the system) and $m_{i>0} = m_i \sigma_{i-1} / \sigma_i$. Noting that $\mathbf{p}_0$ 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 $\mathbf{r}_1 = \mathbf{R}_1 - \mathbf{R}_0$. On the other hand, the two terms involving the planet are expanded in powers of $\eta$. Collecting the terms,

$$H_0 = \frac{\mathbf{p}_1^2}{2m_1} + \frac{\mathbf{p}_2^2}{2m_2} - \frac{Gm_1\sigma_1}{r_1} - \frac{Gm_2\sigma_2}{r_2}, \qquad (B.7)$$

$$H_1 = -\frac{1}{2} \frac{\sigma_2}{\sigma_1} \frac{Gm_1 m_2}{r_2} \left[ \frac{3}{2} \left( \frac{\mathbf{r}_1 \cdot \mathbf{r}_2}{r_2^2} \right)^2 - \left( \frac{\mathbf{r}_1}{r_2} \right)^2 \right] + \bigcirc(\eta^3). \quad (B.8)$$

One can see that $H_1$ is of order $\eta^2$ compared to the last term of Eq. (B.7) and of order $\eta\epsilon$ compared to the penultimate term, which justifies it being treated as a perturbation. We also note that $\sigma_2/\sigma_1 = 1 + \bigcirc(\epsilon)$, 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_0$ is recognised as the Hamiltonian of a mass $m_1$ orbiting around a central mass $\sigma_1$ and a mass $m_2$ independently orbiting around a central mass $\sigma_2$. Thus the solution is exactly given by two Keplerian orbits.

## Appendix B.1: Timing model corresponding to the dominant Keplerian Hamiltonian

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

$$\Delta_R = -\frac{\mathbf{n}_\odot \cdot \mathbf{R}_0}{c}, \qquad (B.9)$$

where $\mathbf{n}_\odot$ 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 $\mathbf{R}_0$. A complete treatment could formally be obtained by replacing $\mathbf{R}_0$ by $\mathbf{R}_p = \mathbf{R}_0 + \delta\mathbf{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,

$$\mathbf{R}_0 = -\frac{M_1}{\sigma_1} \mathbf{r}_1 - \frac{M_2}{\sigma_2} \mathbf{r}_2, \qquad (B.10)$$

where according to the exact Hamiltonian $H_0 \mathbf{r}_{1,2}$ follow Keplerian orbits of the form

$$\mathbf{r}_i = a_i \left( \frac{\cos E_i - e_i}{\sqrt{1 - e_i^2} \sin E_i} \right)_{F_i}, \qquad (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 eccentricity anomaly and $e_i$ the eccentricity. Kepler's equation relates $E_i$ to the mean anomaly $M_i = 2\pi(t - T_i)/P_i$, where $P_i$ is the orbital period and $T_i$ the time of periastron passage,

$$E_i - e_i \sin E_i = M_i. \qquad (B.12)$$

In order to express the Keplerian Rømer delay one needs to express both $\mathbf{n}_\odot$ and $\mathbf{R}_0$ in the same frame. We choose a frame where $\mathbf{n}_\odot$ is the third basis vector. We call it the observer's frame

$F_0$ (who would be standing at the Solar System barycentre). Orbital inclinations $i_j$ are then defined relative to $\mathbf{n}_\odot$ and longitudes of ascending nodes $\Omega_j$ define a rotation around $\mathbf{n}_\odot$. Together with the arguments of periastron $\omega_j$ these angles define the rotations relating the observer's frame and $F_j$,

$$\mathbf{r}_j|_{F_0} = R_{\Omega_j} R_{i_j} R_{\omega_j} \mathbf{r}_j|_{F_j}, \qquad (B.13)$$

where $R_{\Omega_j}$, $R_{i_j}$, $R_{\omega_j}$ are rotation matrices (see e.g. Beutler (2004)), and $\mathbf{r}_j|_F$ refers to the coordinate of $\mathbf{r}_j$ in frame F. From there, the formula expressing the Keplerian Rømer delay $\Delta_j = -\mathbf{n}_\odot \cdot \mathbf{r}_j/c$ is well known (e.g. Lyne & Graham-Smith 2012; Hobbs et al. 2006),

$$\Delta_j = -\frac{a_j \sin i_j}{c} \left( \sqrt{1 - e_j^2} \sin E_j \cos \omega_j + \cos E_j \sin \omega_j - e_j \sin \omega_j \right), \qquad (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 $\Omega_j$ and that $a_j$ cannot be separated from $\sin i_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

$$\Delta_R = - \sum_{j=1}^{2} \frac{M_j}{\sigma_j} \Delta_j. \qquad (B.15)$$

## Appendix B.2: Perturbation

We now turn to the perturbative treatment of $H_1$, 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

$$\Delta_j = \Delta_j^{(0)} + \delta\Delta_j, \qquad (B.16)$$

where $\Delta_j^{(0)}$ is given by Eq. (B.14), and

$$\delta\Delta_j = \delta \left( -\frac{a_j \sin i_j}{c} \left( \sin(n_j t + \phi_j) \right) \right), \qquad (B.17)$$

which denotes the variation of the leading order term in eccentricity of Eq. (B.14). For convenience we introduced $\phi_j = \omega_j - n_j T_j$ such that at leading order $n_j t + \phi_j = E_j + \omega_j$.

The outer-binary term in Eq. (B.15) is dominant since $r_2 M_2 / \sigma_2 \sim (\epsilon/\eta) r_1 M_1 / \sigma_1$ and in the present case $\epsilon \ll \eta$. Therefore we focus on the perturbation of $\Delta_1$. Without loss of generality we choose a frame such that $\Omega_2 = 0$ such that $\Omega_1$ is defined relative to the longitude of ascending node of the planet.

In order to get rid of periodic terms we average $H_1$ over $M_1$, being understood that $P_1 \gg P_2$, and obtain $\bar{H}_1 = (2\pi)^{-1} \int_0^{2\pi} H_1 dM_1$. Although straightforward, this operation is very lengthy. We used the formal calculus software Sage (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 $\omega_{1,2}$, $T_1$, $E_1$, $e_{1,2}$, $\Omega_2$. The relevant Laplace-Lagrange variational equations are

$$ \dot{i}_1 = -\frac{1}{n_1 a_1^2 \sin i_1} \frac{\partial H_1}{\partial \Omega_1}, \qquad (B.18) $$

$$ \dot{\Omega}_1 = \frac{1}{n_1 a_1^2 \sin i_1} \frac{\partial H_1}{\partial i_1}, \qquad (B.19) $$

$$ \dot{\phi}_1 = -\frac{1}{n_1 a_1} \left( \frac{1}{a_1 \tan i_1} \frac{\partial H_1}{\partial i_1} + 2 \frac{\partial H_1}{\partial a_1} \right), \qquad (B.20) $$

$$ \dot{a}_1 = \dot{n}_1 = \dot{e}_1 = 0. \qquad (B.21) $$

We seek solutions of the form $x(t) = x^{(0)} + \delta x(t)$ where $x^{(0)} \in \{a_1, n_1, i_1, \Omega_1, \phi_1\}$ are the constant orbital elements that are solution of the unperturbed system. To leading order in perturbation we may thus substitute $x \to x^{(0)}$ in the right-hand side of Eqs.(B.18)-(B.21) and integrate in order to obtain $\delta x = \int \dot{x} dt$. In order to focus on the main contributions we keep only terms diverging with time, that is terms $\propto E_2$. In the following, we drop the (0) exponent for readability since orbital elements always refer to the unperturbed solution, that is $x = x^{(0)}$ unless otherwise stated.

Substituting $x \to x + \delta x(t)$ into Eq. (B.17) we obtained an expression of the form

$$ \delta\Delta_1 = \alpha t (\gamma_c \cos E_1 + \gamma_s \sin E_1), \qquad (B.22) $$

where $E_1 = n_1(t - T_1)$ at leading order in eccentricity and

$$ \alpha = \frac{3 G m_2 a_1}{8 c n_1 a_2^3}, \qquad (B.23) $$

$$ \gamma_s = -\sin \Omega_1 \cos i_1 \sin i_2 \qquad (\cos i_1 \cos i_2 + \cos \Omega_1 \sin i_1 \sin i_2), \qquad (B.24) $$

$$ \gamma_c = \sin i_1 \left( \cos^2 i_1 \cos^2 i_2 - \frac{10}{3} \right) \qquad (B.25) \\ + \cos \Omega_1 \sin i_2 \cos i_2 \cos i_1 (1 + 2 \sin^2 i_1) \\ + (\cos \Omega_1 \sin i_2)^2 \sin i_1 (1 + \sin^2 i_1). $$

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 $\Delta_1$, 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<fnref n="11" />. The result is shown in Fig.1, where the Keplerian component $\Delta_1^{(0)}$ 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 $\alpha\gamma_c$, $\alpha\gamma_s$ can in principle be fitted independently due to the Fourier decomposition theorem. Since $\alpha$ only occurs as a common scaling in Eq.(B.22), we focus on $\gamma_c$, $\gamma_s$. If their relation to $\Omega_1$, $i_2$ 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 $\gamma_c(\Omega_1, i_2)$, $\gamma_s(\Omega_1, i_2)$ for $i_1 = 39$ deg, which approximately corresponds to the observed inclination of the outer binary. We see that although not degenerate, $\gamma_c$ is about five times larger than $\gamma_s$ for this particular inclination $i_1$. Given the overall weakness of the signal, it might therefore be difficult to constrain accurately both $i_1$ and $\omega_2$ in this particular case.

<figure><img src="image_13.png" /><figcaption>Fig. B.1. Maps representing the values of $\gamma_c$ (left) and $\gamma_s$ (right) as a function of $(\Omega_1, i_2)$ for $i_1 = 39$ deg.</figcaption></figure>

## 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.

<footnote n="11">¹¹ Implemented in the Scipy least_square routine.</footnote>
<table><caption>**Table C.1.** All fits carried out for each model assuming GR.</caption><thead><tr><th>Model</th><th>Npar</th><th>χ<sup>2</sup></th><th>R<sub>χ</sub><sup>2</sup></th><th>ΔAIC</th><th>ΔBIC</th></tr></thead><tbody><tr><td>PL2/Kepl1</td><td>30</td><td></td><td></td><td></td><td></td></tr><tr><td>mean</td><td></td><td>15531.9</td><td>1.24814</td><td>15591.9</td><td>15814.8</td></tr><tr><td>max</td><td></td><td>8.09</td><td>0.00065</td><td>8.09</td><td>8.09</td></tr><tr><td>refit max</td><td></td><td>4.64</td><td>0.000373</td><td>4.64</td><td>4.64</td></tr><tr><td>refit mean</td><td></td><td>-0.915</td><td>-7.35e-05</td><td>-0.915</td><td>-0.915</td></tr><tr><td>PL3</td><td>31</td><td></td><td></td><td></td><td></td></tr><tr><td>mean</td><td></td><td>15543.4</td><td>1.24917</td><td>15605.4</td><td>15835.8</td></tr><tr><td>max</td><td></td><td>7.93</td><td>0.000637</td><td>7.93</td><td>7.93</td></tr><tr><td>refit max</td><td></td><td>-31.5</td><td>-0.00253</td><td>-31.5</td><td>-31.5</td></tr><tr><td>refit mean</td><td></td><td>-30.7</td><td>-0.00247</td><td>-30.7</td><td>-30.7</td></tr><tr><td>PL3DM10</td><td>39</td><td></td><td></td><td></td><td></td></tr><tr><td>mean</td><td></td><td>15496.8</td><td>1.24622</td><td>15574.8</td><td>15864.6</td></tr><tr><td>max</td><td></td><td>-0.509</td><td>-4.09e-05</td><td>-0.509</td><td>-0.509</td></tr><tr><td>refit max</td><td></td><td>-1.87</td><td>-0.000151</td><td>-1.87</td><td>-1.87</td></tr><tr><td>refit mean</td><td></td><td>-2.22</td><td>-0.000179</td><td>-2.22</td><td>-2.22</td></tr><tr><td>PL4</td><td>32</td><td></td><td></td><td></td><td></td></tr><tr><td>mean</td><td></td><td>15569.8</td><td>1.25139</td><td>15633.8</td><td>15871.6</td></tr><tr><td>max</td><td></td><td>7.56</td><td>0.000608</td><td>7.56</td><td>7.56</td></tr><tr><td>refit max</td><td></td><td>-50.3</td><td>-0.00404</td><td>-50.3</td><td>-50.3</td></tr><tr><td>refit mean</td><td></td><td>-53.2</td><td>-0.00427</td><td>-53.2</td><td>-53.2</td></tr><tr><td>PL4 bis</td><td>32</td><td></td><td></td><td></td><td></td></tr><tr><td>mean</td><td></td><td>15558.3</td><td>1.25047</td><td>15622.3</td><td>15860.1</td></tr><tr><td>max</td><td></td><td>9.62</td><td>0.000773</td><td>9.62</td><td>9.62</td></tr><tr><td>refit mean</td><td></td><td>-38.2</td><td>-0.00307</td><td>-38.2</td><td>-38.2</td></tr><tr><td>refit max</td><td></td><td>-36</td><td>-0.00289</td><td>-36</td><td>-36</td></tr><tr><td>PL5</td><td>33</td><td></td><td></td><td></td><td></td></tr><tr><td>mean</td><td></td><td>15535.8</td><td>1.24876</td><td>15601.8</td><td>15847.1</td></tr><tr><td>max</td><td></td><td>4.82</td><td>0.000388</td><td>4.82</td><td>4.82</td></tr><tr><td>refit max</td><td></td><td>-24.2</td><td>-0.00195</td><td>-24.2</td><td>-24.2</td></tr><tr><td>refit mean</td><td></td><td>-29.3</td><td>-0.00236</td><td>-29.3</td><td>-29.3</td></tr><tr><td>planet</td><td>32</td><td></td><td></td><td></td><td></td></tr><tr><td>max</td><td></td><td>15520.9</td><td>1.24746</td><td>15584.9</td><td>15822.7</td></tr><tr><td>mean</td><td></td><td>520</td><td>0.0418</td><td>520</td><td>520</td></tr><tr><td>refit mean</td><td></td><td>-1.76</td><td>-0.000141</td><td>-1.76</td><td>-1.76</td></tr><tr><td>refit max</td><td></td><td>-6.56</td><td>-0.000527</td><td>-6.56</td><td>-6.56</td></tr></tbody></table>

**Notes.** The reduced $\chi^2$, $R_{\chi}^2$, is equal to $\chi^2$/Ndof where Ndof = 13534 - Npar. AIC and BIC are the Akaike and Bayesian information criteria, respectively. For each model the first line is the reference fit subtracted from all the other lines. The last line is the best fit of each given model that appears in Table 2 (except for 'PL4 bis'). The label 'mean' indicates that it is the mean of the MCMC sample, 'max' shows that it is the element in the sample with the highest likelihood, and 'refit' means that a local fit using the *Minuit* library (see sec.4.1) from either the 'mean' or 'max' solution has been performed. For model 'PL4', 'bis' refers to a separate MCMC run (see text).

## 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 $\tilde{f}, \tilde{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.
<figure><img src="image_14.png" /><figcaption>**Fig. D.1.** Correlation plot of the posterior distribution function of the PL3 model using the corresponding MCMC-generated sample, and restricted to *f*, *f* as well as PL3 specific parameters.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. D.2.** Correlation plot of the posterior distribution function of the PlanetGR model using the corresponding MCMC-generated sample, and restricted to $\tilde{f}, \tilde{f}$ as well as Planet specific parameters.</figcaption></figure>
Appendix E: Additional material

<figure><img src="image_16.png" /><figcaption>**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\text{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/P_O \approx 1.5 \times 10^{-2}\text{day}^{-1}$ where $P_O$ is the outer binary period, we show the periodograms of the planet (dashed orange), Kepler (dash-dotted red), PL2/Kep11 (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.</figcaption></figure>
<table><caption>**Table E.1.** Mean values of the MCMC runs of the PL3 and Planet models with their 68% median confidence intervals.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>PL3GR</th><th>PlanetGR</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center;">Fixed values</td></tr><tr><td>Reference epoch (MJD)</td><td>T<sub>ref</sub></td><td>56492</td><td></td></tr><tr><td>Position epoch (MJD)</td><td>T<sub>pos</sub></td><td>57205</td><td></td></tr><tr><td colspan="4" style="text-align: center;">Fitted values</td></tr><tr><td>Right ascension</td><td>α</td><td>3<sup>h</sup>37<sup>m</sup>43<sup>s</sup>.82700(43)<sup>+66</sup><sub>-67</sub></td><td>3<sup>h</sup>37<sup>m</sup>43<sup>s</sup>.82700(94)<sup>+63</sup><sub>-64</sub></td></tr><tr><td>Declination</td><td>δ</td><td>17<sup>deg</sup>15<sup>m</sup>14<sup>s</sup>.817(63)<sup>+38</sup><sub>-37</sub></td><td>17<sup>deg</sup>15<sup>m</sup>14<sup>s</sup>.817(75)<sup>+40</sup><sub>-33</sub></td></tr><tr><td>Distance (kpc)</td><td>d</td><td>1160<sup>+53</sup><sub>-51</sub></td><td>1182<sup>+53</sup><sub>-55</sub></td></tr><tr><td>Right-ascension proper motion (mas yr<sup>-1</sup>)</td><td>μ<sub>α</sub></td><td>6.1(37)<sup>+70</sup><sub>-68</sub></td><td>5.9(82)<sup>+70</sup><sub>-80</sub></td></tr><tr><td>Declination proper motion (mas yr<sup>-1</sup>)</td><td>μ<sub>δ</sub></td><td>-5. (34)<sup>+25</sup><sub>-26</sub></td><td>-4. (83)<sup>+29</sup><sub>-25</sub></td></tr><tr><td>Radial proper motion (mas yr<sup>-1</sup>)</td><td>μ<sub>d</sub></td><td>5. (63)<sup>+31</sup><sub>-32</sub></td><td>5. (52)<sup>+32</sup><sub>-28</sub></td></tr><tr><td>Dispersion measure (pc cm<sup>-3</sup>)</td><td>DM</td><td>21.316(46)<sup>+12</sup><sub>-12</sub></td><td>21.316(46)<sup>+11</sup><sub>-11</sub></td></tr><tr><td>Dispersion measure variation (pc cm<sup>-3</sup> yr<sup>-1</sup>)</td><td>DM'</td><td>(1.7)<sup>+2.1</sup><sub>-2.1</sub> × 10<sup>-5</sup></td><td>(1.6)<sup>+1.9</sup><sub>-2.0</sub> × 10<sup>-5</sup></td></tr><tr><td>Re-scaled spin frequency (Hz)</td><td>f&#773;</td><td>365.953337902243(85)<sup>+81</sup><sub>-81</sub></td><td>365.9533379022(56)<sup>+11</sup><sub>-11</sub></td></tr><tr><td>Re-scaled spin frequency derivative (10<sup>-15</sup> Hz s<sup>-1</sup>)</td><td>f&#773;'</td><td>-2.35631(40)<sup>+48</sup><sub>-45</sub> × 10<sup>-15</sup></td><td>-2.356(45)<sup>+11</sup><sub>-12</sub> × 10<sup>-15</sup></td></tr><tr><td colspan="4" style="text-align: center;">Inner orbit I of pulsar p w.r.t. inner-binary centre of mass b</td></tr><tr><td>Orbital period (d)</td><td>P<sub>1</sub></td><td>1.62939900(75)<sup>+42</sup><sub>-42</sub></td><td>1.62939900(83)<sup>+38</sup><sub>-41</sub></td></tr><tr><td>Projected semi-major axis (lt-s)</td><td>a<sub>p</sub> sin i<sub>1</sub></td><td>1.2175276(17)<sup>+56</sup><sub>-56</sub></td><td>1.2175276(25)<sup>+54</sup><sub>-56</sub></td></tr><tr><td>Inclination offset (°)</td><td>δi = i<sub>1</sub> - i<sub>0</sub></td><td>1.2<sup>+2.3</sup><sub>-2.4</sub> × 10<sup>-3</sup></td><td>8<sup>+21</sup><sub>-23</sub> × 10<sup>-4</sup></td></tr><tr><td>Laplace-Lagrange</td><td>e<sub>1</sub> cos ω<sub>p</sub></td><td>6.937(38)<sup>+40</sup><sub>-40</sub> × 10<sup>-4</sup></td><td>6.937(32)<sup>+37</sup><sub>-42</sub> × 10<sup>-4</sup></td></tr><tr><td>Laplace-Lagrange</td><td>e<sub>1</sub> sin ω<sub>p</sub></td><td>-8.59(37)<sup>+40</sup><sub>-41</sub> × 10<sup>-5</sup></td><td>-8.59(42)<sup>+36</sup><sub>-35</sub> × 10<sup>-5</sup></td></tr><tr><td>Time of ascending node (MJD)</td><td>t<sub>ascp</sub></td><td>55917.15894(93)<sup>+19</sup><sub>-19</sub></td><td>55917.15894(89)<sup>+19</sup><sub>-17</sub></td></tr><tr><td>Long. of asc. nodes offset (°)</td><td>δΩ = Ω<sub>b</sub> - Ω<sub>p</sub></td><td>(1.7)<sup>+1.3</sup><sub>-1.3</sub> × 10<sup>-4</sup></td><td>(1.5)<sup>+1.2</sup><sub>-1.2</sub> × 10<sup>-4</sup></td></tr><tr><td colspan="4" style="text-align: center;">Outer orbit O of b w.r.t. triple-system centre of mass t</td></tr><tr><td>Orbital period (d)</td><td>P<sub>O</sub></td><td>327.25512(55)<sup>+10</sup><sub>-10</sub></td><td>327.25512(61)<sup>+10</sup><sub>-11</sub></td></tr><tr><td>Projected semi-major axis (lt-s)</td><td>a<sub>b</sub> sin i<sub>O</sub></td><td>74.672344(35)<sup>+22</sup><sub>-23</sub></td><td>74.672344(31)<sup>+21</sup><sub>-22</sub></td></tr><tr><td>Co-projected semi-major axis (lt-s)</td><td>a<sub>b</sub> cos i<sub>O</sub></td><td>91.4(10)<sup>+37</sup><sub>-37</sub></td><td>91.4(02)<sup>+35</sup><sub>-32</sub></td></tr><tr><td>Laplace-Lagrange</td><td>e<sub>O</sub> cos ω<sub>b</sub></td><td>0.03511469(40)<sup>+28</sup><sub>-27</sub></td><td>0.03511469(43)<sup>+28</sup><sub>-27</sub></td></tr><tr><td>Laplace-Lagrange</td><td>e<sub>O</sub> sin ω<sub>b</sub></td><td>-0.00352492(40)<sup>+19</sup><sub>-18</sub></td><td>-0.00352492(22)<sup>+19</sup><sub>-18</sub></td></tr><tr><td>Time of ascending node (MJD)</td><td>t<sub>ascb</sub></td><td>56230.195329(47)<sup>+83</sup><sub>-82</sub></td><td>56230.195329(01)<sup>+79</sup><sub>-78</sub></td></tr><tr><td>Longitude of ascending node (°)</td><td>Ω<sub>b</sub></td><td>145.9<sup>+1.4</sup><sub>-1.5</sub></td><td>148.0<sup>+1.7</sup><sub>-1.5</sub></td></tr><tr><td>Inner mass ratio</td><td>m<sub>i</sub>/m<sub>p</sub></td><td>0.1373(76)<sup>+12</sup><sub>-12</sub></td><td>0.1373(77)<sup>+12</sup><sub>-12</sub></td></tr><tr><td>ToA uncertainty re-scaling</td><td>EF<sub>AC</sub></td><td>1.11(77)<sup>+76</sup><sub>-73</sub></td><td>1.11(72)<sup>+75</sup><sub>-79</sub></td></tr><tr><td colspan="4" style="text-align: center;">Derived values</td></tr><tr><td>Parallel proper motion (km s<sup>-1</sup>)</td><td>V<sub>||</sub></td><td>29.77<sup>+0.90</sup><sub>-1.1</sub></td><td>29. (72)<sup>+76</sup><sub>-83</sub></td></tr><tr><td>Plane-of-sky proper motion (km s<sup>-1</sup>)</td><td>V<sub>⊥</sub></td><td>43.0<sup>+2.2</sup><sub>-2.0</sub></td><td>41.4<sup>+2.1</sup><sub>-2.2</sub></td></tr><tr><td>Spin frequency (Hz)</td><td>f</td><td>365.9533629(82)<sup>+13</sup><sub>-13</sub></td><td>365.9533629(80)<sup>+12</sup><sub>-12</sub></td></tr><tr><td>Spin frequency derivative (10<sup>-15</sup> Hz s<sup>-1</sup>)</td><td>f&#773;</td><td>-2.29(06)<sup>+43</sup><sub>-40</sub> × 10<sup>-15</sup></td><td>-2.29(67)<sup>+41</sup><sub>-43</sub> × 10<sup>-15</sup></td></tr><tr><td colspan="4" style="text-align: center;">Inner orbit I of p/b</td></tr><tr><td>Semi-major axis (lt-s)</td><td>a<sub>p</sub></td><td>1.924(46)<sup>+48</sup><sub>-48</sub></td><td>1.924(39)<sup>+46</sup><sub>-43</sub></td></tr><tr><td>Orbital inclination (°)</td><td>i<sub>1</sub></td><td>39.2(47)<sup>+11</sup><sub>-11</sub></td><td>39.2(48)<sup>+10</sup><sub>-11</sub></td></tr><tr><td>Orbital eccentricity</td><td>e<sub>1</sub></td><td>6.990(40)<sup>+40</sup><sub>-40</sub> × 10<sup>-4</sup></td><td>6.990(36)<sup>+37</sup><sub>-42</sub> × 10<sup>-4</sup></td></tr><tr><td>Longitude of periastron (°)</td><td>ω<sub>p</sub></td><td>-7.06(16)<sup>+33</sup><sub>-33</sub></td><td>-7.06(20)<sup>+29</sup><sub>-28</sub></td></tr><tr><td>Time of periastron passage (MJD)</td><td>t<sub>p1</sub></td><td>55917.1269(88)<sup>+15</sup><sub>-15</sub></td><td>55917.1269(86)<sup>+13</sup><sub>-13</sub></td></tr><tr><td>Longitude of asc. node (°)</td><td>Ω<sub>p</sub></td><td>145.9<sup>+1.4</sup><sub>-1.5</sub></td><td>148.0<sup>+1.7</sup><sub>-1.5</sub></td></tr><tr><td colspan="4" style="text-align: center;">Outer orbit O of b w.r.t. triple-system centre of mass t</td></tr><tr><td>Semi-major axis (lt-s)</td><td>a<sub>b</sub></td><td>1.180(33)<sup>+28</sup><sub>-29</sub> × 10<sup>2</sup></td><td>1.180(27)<sup>+27</sup><sub>-25</sub> × 10<sup>2</sup></td></tr><tr><td>Orbital inclination (°)</td><td>i<sub>O</sub></td><td>39.2(45)<sup>+11</sup><sub>-11</sub></td><td>39.2(48)<sup>+10</sup><sub>-10</sub></td></tr><tr><td>Orbital eccentricity</td><td>e<sub>O</sub></td><td>0.03529117(20)<sup>+29</sup><sub>-29</sub></td><td>0.03529117(21)<sup>+29</sup><sub>-28</sub></td></tr><tr><td>Longitude of periastron (°)</td><td>ω<sub>b</sub></td><td>-5.73232(85)<sup>+27</sup><sub>-27</sub></td><td>-5.73232(56)<sup>+27</sup><sub>-27</sub></td></tr><tr><td>Time of periastron passage (MJD)</td><td>t<sub>pO</sub></td><td>56224.98440(20)<sup>+19</sup><sub>-19</sub></td><td>56224.98440(42)<sup>+19</sup><sub>-19</sub></td></tr><tr><td>Pulsar mass (M<sub>☉</sub>)</td><td>m<sub>p</sub></td><td>1.43(82)<sup>+12</sup><sub>-12</sub></td><td>1.43(80)<sup>+12</sup><sub>-11</sub></td></tr><tr><td>Inner-companion mass (M<sub>☉</sub>)</td><td>m<sub>i</sub></td><td>0.197(58)<sup>+15</sup><sub>-15</sub></td><td>0.197(55)<sup>+15</sup><sub>-14</sub></td></tr><tr><td>Outer-companion mass (M<sub>☉</sub>)</td><td>m<sub>o</sub></td><td>0.410(18)<sup>+31</sup><sub>-31</sub></td><td>0.410(12)<sup>+31</sup><sub>-29</sub></td></tr></tbody></table>

**Notes.** Error bars apply to the digits between parenthesis and delimit the 68% median confidence region, namely the interval excluding 16% of the distribution above and below it. The central value quoted is the best-fit value. Upper-case indices I, O refer to the inner and outer binaries respectively. Lower-case indices p, i, o refer to the pulsar, inner white dwarf, and outer white dwarf, respectively; b refers to the inner-binary centre of mass. Data spans MJD 56492 – 59480. Solar System ephemeris is DE430. Model-specific parameters are reported separately in Table 1.