© ESO 2020

1. Introduction

Among the fundamental interactions of nature, gravity is unique in attracting all material objects with the same acceleration, at least within current observational precision. This feature of gravity (the universality of free fall, UFF below) was thought by Newton to be a cornerstone of Newtonian mechanics (Newton 1687). Indeed, in the Newtonian theory of gravity, this universal acceleration implies that the inertial mass of a body is always in a fixed proportion to its passive gravitational mass, and is independent of the mass, chemical composition, or the detailed internal structure of the gravitating object. This was presented as an observed physical principle, without a deeper explanation. Newton and many later experimentalists have conducted different experiments to verify UFF, no deviations have been found that are larger than 1.3 × 10⁻¹⁴ (Touboul et al. 2019). This equivalence between the inertial and passive gravitational masses for test particles (defined here as objects with negligible gravitational self-energy) is the so-called weak equivalence principle (WEP).

When thinking about a new theory of gravity that incorporates the laws of special relativity (SR), Einstein had the insight that the gravitational field appears to be absent for a freely falling observer. This was later described by Einstein as the ‘most fortunate thought in my life’ (Renn 2007). This idea, that gravity is equivalent to acceleration, naturally explains the WEP. If the relativity principle applies to this situation, then any observers in a sufficiently small room in a free-falling reference frame are not only unable to determine whether the room is in motion or at rest relative to distant bodies, but they are neither able to determine its rate of acceleration in the gravitational field. This implies that, in the vicinity of the observer the laws of physics are (in very good approximation) given by SR, which means that the Lorentz invariance of SR is obeyed locally (this is the local Lorentz invariance, LLI) and furthermore, that it does not matter where or when an experiment is made (this is known as local position invariance, LPI). The combination of the WEP with LLI and LPI is now known as the Einstein equivalence principle (EEP, Will 2018a). Schiff’s conjecture states that the WEP implies the full EEP for any consistent theory of gravity, for which a strong plausibility argument can be made (see e.g. Will 2018a).

This generalisation of the relativity principle to reference frames in free fall guided Einstein towards general relativity (GR, Einstein 1915). GR and other metric theories of gravity fulfil the EEP in a natural way: in these theories the gravitational attraction is seen as a result of spacetime curvature, which itself originates from the energy, stress and momentum of the masses in a system, determined by the field equations of the theory. This curvature changes the trajectories of test particles moving within the spacetime (their ‘geodesics’) in a unique way that does not depend of the detailed nature of the particles themselves, hence the validity of the WEP. Furthermore, for spatial scales that are small compared to the radius of curvature, the geometry of spacetime necessarily approximates the ‘flat’ Minkowski geometry, hence the LLI and LPI automatically apply to non-gravitational experiments. To rephrase, the EEP is a consequence of a universal coupling between matter and gravity (Damour 2012).

The qualification of ‘non-gravitational’ is key here. If the EEP can be fully extended to gravitational experiments, like the Cavendish experiment, and to objects with large self-gravitational energy, then we have the strong equivalent principle (SEP). This is a crucially important distinction because, while all metric theories of gravity fulfil the EEP, there are suggestive arguments that GR is the only gravity theory in four spacetime dimensions that fully embodies the SEP (Di Casola et al. 2015; Will 2018a)¹.

Therefore, if we are looking for phenomena beyond GR, a promising avenue would be to look for instances of SEP violation. This has an added advantage: if no SEP violation is found, the results of such an experiment can in principle constrain all alternative theories of gravity.

Just as the EEP consists of the WEP, LLI and LPI, the SEP must additionally include gravitational versions of these. Any violations of the LLI and LPI of the gravitational interaction (e.g., the existence of a preferred frame of reference or the location dependence of gravity) have been strongly constrained using pulsar experiments (Shao & Wex 2012, 2013; Shao et al. 2013). In what follows, we focus on the gravitational version of the WEP (GWEP, Will 2018a), which states that the UFF applies not only to test particles, but also to any objects where the gravitational binding energy is important.

For alternative theories of gravity, the gravitational properties of objects generally depend on their amount of self-gravity. This means that at Newtonian level we have a body-dependent effective gravitational constant, G_ab, meaning the acceleration of a body a in the gravitational field of a body b is given by

$$\begin{array}{c}\hfill {\ddot{\mathit{x}}}_a=-G_{\mathit{ab}}{m}_b\frac{{\mathit{r}}_{\mathit{ab}}}{\Vert {\mathit{r}}_{\mathit{ab}}{\Vert }^3}+\mathcal{O}(c^{-2}),\end{array}$$(1)

where m_b denotes the inertial mass of body b, r_ab ≡ x_a − x_b their (coordinate) separation, and c is the speed of light. G_ab depends on the properties of body a and b. In the weak-field limit this can be interpreted, to a good approximation, as a mismatch between the inertial and the gravitational masses of the objects:

$$\begin{array}{c}\hfill G_{\mathit{ab}}={\left(\frac{m_{\mathrm{P}}}{m}\right)}_a{\left(\frac{m_{\mathrm{A}}}{m}\right)}_b{G}_{\mathrm{N}},\end{array}$$(2)

where G_N is the Newtonian gravitational constant, as measured in a Cavendish-type experiment, and m_P and m_A denote the passive and active gravitational mass respectively. For semi-conservative metric theories of gravity that have a conservation of momentum one has only a single gravitational mass m_G ≡ m_P = m_A (Will 2018a). For the remainder of the paper we assume that momentum is conserved in the gravitational interaction, and therefore G_ab = G_ba². More generally, we use the definition G_ab = G_N(1 + Δ_ab) where we denote Δ_ab = Δ_ba as the relative GWEP parameter between two bodies a and b.

If one observes an isolated two-body system without prior knowledge of the masses, then any violation of the GWEP at Newtonian order would be indistinguishable from a re-scaling of the masses due to the symmetry of the equations of motion. This symmetry is broken in presence of a third body. One can then compare the rate at which two self-gravitating objects fall in the field of a third one. This forms the base for a class of GWEP tests that includes Lunar laser ranging (LLR), tests with binary pulsars falling in the gravitational field of our Galaxy, and the test to be discussed in this paper.

In the LLR test, one considers the Earth-Moon system falling in the gravitational field of the Sun. If GWEP is violated, then the Earth, which has a larger fractional gravitational binding energy than the Moon, falls in the Sun’s gravitational field with a slightly different acceleration than the Moon. This causes a polarisation of the Earth-Moon orbit in the direction of the Sun (Nordtvedt 1968). This so called Nortdvedt effect is the gravitational equivalent of the Stark effect, where a strong electric field polarises neutral atoms. It manifests itself as an added small orbital eccentricity vector that precesses in the sky with a period of 1 year, trailing the Sun. The relative Earth-Moon distance can be measured with an accuracy of about 10 cm thanks to the reflectors laid on the Moon by a variety of American and Soviet lunar missions. No Nodtvedt effect has been measured, as predicted by GR, effectively constraining

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{{\mathrm{E}}_{\odot }}-{\mathrm{\Delta }}_{{\mathrm{M}}_{\odot }}\simeq {\left(\frac{m_{\mathrm{G}}}{m}\right)}_{\mathrm{E}}-{\left(\frac{m_{\mathrm{G}}}{m}\right)}_{\mathrm{M}}=(-3.0\pm 5.0)\times {10}^{-14}\end{array}$$(3)

(Hofmann & Müller 2018), which is only about a factor of 10 weaker than the MICROSCOPE limit for WEP, therefore confirming to a high degree that gravitational binding energy falls the same way in an external gravitational field as any other form of energy.

In this test, all the involved bodies are weakly self-gravitating, however, this is especially true for the two ‘proof masses’, the Earth and the Moon: for the Earth ε_grav, E ≡ E_grav, E/m_Ec² = −4.6 × 10⁻¹⁰ (here E_grav, E is the Newtonian gravitational binding energy of the Earth), for the Moon ε_grav, M = −0.2 × 10⁻¹⁰. This means that the LLR experiment only tests the weak-field limit of GWEP. In this limit Eq. (2) implies Δ_ab ≃ Δ_a + Δ_b, where Δ_a ≡ (m_G/m)_a − 1; furthermore, the gravitational binding energy of the bodies relative to their mass is so small that it can only have a very small effect on Δ_a. Within the parametrised post-Newtonian (PPN) formalism for metric theories of gravity

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_a=\eta {\epsilon }_{\mathrm{grav},a},\end{array}$$(4)

where η is the so called Nordtvedt parameter, a combination of several PPN parameters (see Will 2018a, for details). The current limit on the Nordtvedt parameter from LLR is (− 0.2 ± 1.1) × 10⁻⁴.

A violation of GWEP not only affects the dynamics of the Earth-Moon system, but all self-gravitating masses in the Solar System are affected according to Eq. (4). A consequence of this is a shift of the Solar System barycentre (SSB) when modelling planetary ephemerides. Based on data from the MESSENGER mission, Genova et al. (2018) have derived η = (− 6.6 ± 7.2) × 10⁻⁵.

Equation (4) applies to the weak-field limit, that is, the Nordtvedt parameter parametrises GWEP violation to leading order in ε_grav, a ≪ 1. This first order approach is no longer applicable in the strong-field regime of neutron stars. Thus, in the remainder of this article we consider GWEP violations in terms of limits directly on Δ_ab and not on η.

In the Damour-Schäfer test (Damour & Schäfer 1991), one verifies whether the two components of a pulsar – white dwarf system (the first with a very high degree of self-gravity, which allows the detection of strong-field SEP violation) fall with the same acceleration in the field of the Galaxy, which acts as the third body. A violation of the UFF would again cause a polarisation of the orbit of the binary pulsar. At the time of that paper (1991), the timing precision and timing baselines of binary pulsars were relatively small, so the authors proposed a statistical approach to search for this polarising effect in the orbital eccentricities of the known pulsar – white dwarf systems. Following that method, several analyses of the orbital eccentricities have constrained Δ for neutron stars³: Stairs et al. (2005) derive |Δ| < 5.6 × 10⁻³, and Gonzalez et al. (2011) derive |Δ| < 4.6 × 10⁻³ (both being 95% confidence limits). However, the latter limit is derived with the inclusion of a binary pulsar, PSR J1711−4322 that does not fulfil all the necessary criteria for the Damour-Schäfer test (Wex 2014).

This method has several shortcomings, which are listed and discussed in detail by Freire et al. (2012a); the two most important ones are (a) the fact that it cannot detect GWEP violation, only produce statistical upper limits for it and (b) generally, neutron stars with different masses have different values of Δ, this limits the meaningfulness of a general Δ for neutron stars (cf. footnote 25 in Damour 2009).

Apart from the statistical test based on small eccentricities, Damour & Schäfer (1991) have also proposed a test based on a direct measurement of the variation of the orbital eccentricity vector for individual systems, ė, (no matter whether eccentric or not) that results from the polarisation of the binary orbit by the Nordtvedt effect. As discussed by Freire et al. (2012a), this test not only avoids all the shortcomings of the statistical test, but its precision just keeps improving with the precision of the measurement of ė, which improves with time and with better timing instruments. Indeed, they estimated that this test should, for the best timed binaries, yield slightly better Δ values than the statistical test by the mid 2010’s. More recently Zhu et al. (2019) confirmed this by using the ė constraint for the wide orbit of PSR J1713+0747 to derive |Δ| < 2 × 10⁻³ (95% C.L.). Without further assumptions, this limit is strictly speaking only for neutron stars around 1.3 M_⊙, which is the mass of PSR J1713+0747. Recently, this limit has been used to test the UFF of a neutron star towards dark matter (Shao et al. 2018).

Although this test can detect any hypothetical large strong-field deviations of the gravitational properties of neutron stars, the limits on |Δ| are not very constraining because of the weak gravitational field of the Galaxy, which has accelerations of the order of 2 × 10⁻¹⁰ m s⁻² in the solar vicinity. In the case of the LLR test, the polarising gravitational field (that of the Sun) is much stronger (6 × 10⁻³ m s⁻²), however, in that case the Earth and Moon have very small gravitational self energies.

For this reason, Freire et al. (2012a) suggested that the (then) rumoured pulsar in a triple system would combine the best features of both tests. In this experiment, we look for the Nordtvedt effect in an inner binary system consisting of a pulsar and a white dwarf; this system is orbited by a third hierarchical component significantly farther away. As in previous binary pulsar experiments, the pulsar provides the precise tracking and an object with very strong gravitational self energy; the white dwarf provides a test mass with a much smaller gravitational self energy, and finally the third outer component in that system provides a strong (potentially) polarising gravitational field (1.7 × 10⁻³ m s⁻²), as the Sun does for the LLR experiment. The outer component would ideally be a neutron star, as this would provide a qualitatively different test, however, any type of star would already yield a much stronger polarising force than the Galactic gravitational field and therefore either a detection of GWEP violation, or much improved limits on it.

PSR J0337+1715 was discovered in data from the GBT drift-scan survey (Boyles et al. 2013; Lynch et al. 2013). This is a 2.7 ms pulsar in a 1.6 day orbit with a ∼0.2 M_⊙ Helium white dwarf star, this is what we refer to from now on as the inner binary. The outer ∼0.4 M_⊙ white dwarf orbits the inner binary in about 327 days in a low-eccentricity (e = 0.035) orbit. This is the first, and so far the only pulsar confirmed to be in a triple stellar system (Ransom et al. 2014). The two orbits (inner and outer) are nearly coplanar, this and the small observed eccentricities provide important clues for the evolution of the system, which was described in detail by Tauris & van den Heuvel (2014).

The pulsar has very good rotational stability, as usually millisecond pulsars (MSPs) do, and is relatively bright, which allows a very good measurement of the times of arrival of the pulses. This has allowed precise measurements of the varying orbital parameters, and also extremely precise mass measurements for the pulsar and the two white dwarf stars (Ransom et al. 2014). More importantly, the GWEP test was eventually carried out for PSR J0337+1715 by Archibald et al. (2018), yielding |Δ|≤2.6 × 10⁻⁶ (95% C.L.). This represents an improvement of three orders of magnitude over previous pulsar tests and confirmed the power of a MSP in a triple stellar system for testing the GWEP.

This measurement represents a tight constraint on alternative theories of gravity. Archibald et al. (2018) derived constraints on one of the best studied alternatives to GR, the class of mono-scalar-tensor theories described by Damour & Esposito-Farèse (1992; 1993, henceforth DEF gravity). The constraints on the weak-field coupling parameter for these theories (α₀) derived from PSR J0337+1715 significantly improve upon all previous experiments for most of their β₀ space.

The UFF experiment with the PSR J0337+1715 triple system and its results are clearly of great importance. It is, at present, the most powerful test of the GWEP, for either the strong or weak field limits. It is also extremely sensitive to strong-field deviations in the gravitational properties of neutron stars.

This test is of special value because, according to a gravitational analogue of Schiff’s conjecture, it is plausible that the validity of GWEP implies the SEP (Will 2018a); this in turn strongly suggests, according to the arguments mentioned above, that GR is the theory of gravity (see also Will 2014a).

For all these reasons, we find it is important to improve both the precision and reliability of the test. These are the primary objectives of this work. We use fully independent observational data, taken with a wholly different observing system (described in detail in Sect. 2), a completely independent numerical integration of the motion of the system and a different implementation of the determination of the masses and orbital parameters (described in Sect. 3) than those used by Archibald et al. (2018). One of the main differences is, however, that the uncertainties we report are purely statistical; we found no need to postulate the existence of additional systematic effects that can bias Δ. Consequently, this parameter can be self-consistently processed like the others without requiring a special treatment.

The results of our analysis are presented in Sect. 4. Here we discuss not only the parameters we obtain, but also analyse the trends observed in the residuals after subtracting the best model for the system. In Sect. 5, we interpret the Δ constraint, as well as constraints on the post-Newtonian strong-field parameters, within the context of a wide framework of alternative theories of gravity, the Bergmann-Wagoner theories of gravity (Will 2018a). We also derive new constraints on a sub-class of those theories, the Damour-Esposito Farèse (DEF) theories (Damour & Esposito-Farèse 1992, 1993), these new limits are derived in a conservative way that accounts for uncertainties in our knowledge of the equation of state (EOS) for neutron star matter. We finally summarise our findings in Sect. 6.

2. Observations and data reduction

The pulsar J0337+1715 has been regularly observed since July 2013 every 2 or 3 days with the Nançay radio telescope using its L-band receiver at a central frequency of 1484 MHz. The Nançay radio telescope is a meridian Kraus design collector equivalent to a 94 m dish able to conduct ∼1 h 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), an instrument that is able to coherently dedisperse (Hankins & Rickett 1975) a total bandwidth of 512 MHz. It consists of a ROACH1 board (designed by the CASPER group, University of California, Berkeley) providing 128 baseband data streams of 4MHz bandwidth each. The instrument software has many similarities with GUPPI (Green Bank Ultimate Pulsar Processing Instrument, DuPlain et al. 2008) used at the Green Bank Telescope (GBT). A cluster of four nodes hosting eight GTX280/285 Graphics Processing Units (GPUs) is used to coherently dedisperse and fold the data in real-time.

The real-time folding process uses a pulsar period coming from a simple model with two non-interacting Keplerian orbits over short 15 s sub-integrations. A single standard timing parameter file in TEMPO format⁴ containing this pulsar timing model is used for all the observations. The full frequency and time resolution daily pulsar profiles are stored in PSRFITS archives (Hotan et al. 2004)⁵. A ∼3 Hz pulsed noise diode is fired for ∼10 s at the start of each observation to conduct a simple calibration accounting for gain and phase differences between the two polarisations, as implemented in the SINGLEAXIS polarisation calibration of PSRCHIVE.

As the pulsar period model used to fold in real time is not strictly satisfactory, it is necessary to properly phase-shift all the archived individual profiles. A posteriori, the drifts within individual sub-integrations were statistically smaller than the mean ToA uncertainty and characterised by an rms of 0.78 μs, thus validating the parameters of our sub-integrations. In an iterative way, measured times of arrival are used to derive a pulsar timing model which is used to improve the times of arrival and so on. Practically, the numerically derived pulsar timing model is used to provide ‘theoretical’ barycentric arrival times for each of the 15 s sub-integrations. A code transforms those barycentric arrival times into a simple daily TEMPO parameter file with rotation rate described by only a frequency and its first three derivatives around an epoch corresponding to the middle of the observation (F0, F1, F2, F3 and PEPOCH). This polynomial predicts the barycentric rotational phases within 5 ns at maximum. These daily parameter files are then used to re-align the pulse profiles within their corresponding archives. A profile ‘template’, built with more than one thousand observations (see Fig. 1, top), is used for determining the topocentric pulse times of arrival (ToAs) in the following way. After integrating profiles over 128 MHz and 20 min, the times of arrival are estimated using the pat function from PSRCHIVE within the Fourier domain with Markov chain Monte Carlo (FDM) method. The total bandwith of 512 MHz was thus split in four sub-bands in order to be able to fit for variations in the dispersion measure (DM) representing the integrated electron column density along the line-of-sight during the numerical fit. The ToA uncertainties as reported by pat are characterised by a mean of 2.15 μs and a median of 1.89 μs. A pulse profile typical of a good observation, characterised by an uncertainty of 1.15 μs, is shown in Fig. 1 (bottom). The goodness of fit as reported by pat can give a sense of the differences between the template profile and the profiles used to derive ToAs. The goodness values are characterised by a mean of 1.05 (with a median of 1.04) and an rms of 0.12 with 99% of the values between ±3σ (0.69 to 1.41). The rather low signal to noise ratio of the PSR J0337+1715 profiles observed at Nançay prevents the detection of subtle effects of incorrect polarisation calibration on the ToA determination. In this work, we use a dataset (see footnote 7) of 9303 ToAs divided in four 128 MHz bands observed between MJD 56492 and MJD 58761 (July 2013 and October 2019).

Fig. 1. Top: template pulse profile of PSR J0337+1715 used for timing. 450 h of observations conducted with the Nançay Radio Telescope were integrated over the frequency range 1230–1742 MHz. Bottom: pulse profile obtained on October 4, 2014. When profiles are ToA uncertainty ranked, this one is at the 10th percentile of the lowest uncertainties, typical of a good observation.

3. NUmerical TIming MOdel: NUTIMO

For the description of the orbits of binary pulsars, timing programmes such as the aforementioned TEMPO and also TEMPO2 (Edwards et al. 2006; Hobbs et al. 2006) rely on existing analytical models to calculate the times of arrival with nanosecond accuracy (such as, e.g. the DD and DDGR models, Damour & Deruelle 1986). These models are built from precise analytical solutions of the equations of motion (for the examples above these were derived by Damour & Deruelle 1985). However, no such solution is yet available for a triple system where the three masses are of comparable size and experience moderately strong mutual interactions. Therefore, and similarly to Ransom et al. (2014) and Archibald et al. (2018), we perform a high precision numerical integration of the equations of motion, which we subsequently use to calculate the delays.

The equations of motion we use are accurate up to first post-Newtonian order (1PN), that is, include the first-order terms of an expansion of GR in the small parameter $ϵ\sim \frac{\mathit{Gm}}{a{c}^2}\sim \frac{{\mathit{v}}^2}{c^2}\lesssim 5\times {10}^{-7}$ where m, v and a are characteristic mass, velocity and length scales of the system and c the speed of light. These corrections are absolutely necessary because they translate into a relative acceleration Δa/a ∼ ϵ which is of similar magnitude as a potential SEP violation (see above). In addition, 1PN corrections are responsible for effects that accumulate over time such as the well-known relativistic precession of periastron. On the other hand, second order corrections can be safely ignored since the same line of reasoning predicts that even a cumulative effect such as gravitational wave radiation cannot account for more than a nanosecond within the current span of our observations.

We use the 1PN generic strong-field framework of Will (1993) and Damour & Taylor (1992) which parametrises almost the entire class of ‘fully conservative’ Lagrangian-based theories of gravity (without preferred location effects⁶) based on a modified Einstein-Infeld-Hoffmann approach (see details in Appendix A). In this framework, in the most generic case, one has three parameters Δ_ab at the Newtonian and 12 strong-field parameters at the post-Newtonian level. All these parameters depend on the structure of the individual bodies. The 12 post-Newtonian parameters generalise the parametrised post-Newtonian (PPN) β_PPN and γ_PPN (Will 1993) to the regime of strongly self gravitating masses.

In the PSR J0337+1715 system we have only one strongly self-gravitating body with ε_grav ∼ 0.1, the pulsar, while the two white dwarfs are weak-field objects with ε_grav ≲ 10⁻⁴. That generally leads to a significant reduction of the number of strong-field parameters relevant for the orbital dynamics of the PSR J0337+1715 system. In fact, on the Newtonian level there is only one Δ_ab, which we simply denote by Δ. Among the post-Newtonian terms, as we discuss in detail within a theory based framework in Sect. 5, there remain three strong field parameters which are a priori unconstrained by Solar System experiments and limits on Δ already imposed by the ‘Newtonian’ dynamics: ${\overline{\beta }}_{\mathrm{p}}$, ${\overline{\beta }}_{\mathrm{pp}}$, ${\overline{\beta }}^{\mathrm{p}}$. Since the limits we find for these parameters in Sect. 4 are many orders of magnitude weaker than limits inferred indirectly from binary pulsar experiments, at least within our theory based framework, we primarily consider a model where the 1PN strong-field parameters are set to their GR values that is, zero. This practically corresponds to using priors arising from a combination of Solar System and binary pulsar limits at the post-Newtonian level when estimating Δ.

Our specially developed software, NUTIMO⁷, solves numerically the 3-body equations of motion at 1PN (see Appendix A.1 and particularly Eq. (A.2)) before computing propagation and relativistic delays. All the geometrical delays are taken into account up to first order in L/d where d is the distance to the system and L ≪ d is any other length scale of the problem. In other words, the code computes the so-called Rœmer, Kopeikin and Shklovskii delays (Shklovskii 1970; Kopeikin 1996), and adds an extra second order correction, that is, L²/d², for the latter (the only second order correction that may become important with time, see e.g. Voisin 2017). We note that Kopeikin and Shklovskii delays were not included in previous works (Archibald et al. 2018; Ransom et al. 2014). The former allows us to measure the longitude of ascending node of the outer orbit and might be important because it accounts for systematic effects at the Earth orbital frequency which is close to the outer orbital frequency. We do not expect the latter to significantly affect the results of the fit but it allows us to derive the intrinsic pulsar spin parameters which would otherwise absorb this effect (see below). We caution that the intrinsic spin parameters we report in Table 1 are still biased by the effect of Galactic acceleration which amounts to approximately 25% of the Shklovskii correction (according to the Galactic model of McMillan 2017). Relativistic delays include time dilation between the frame of the pulsar and the frame of the observer, namely the so-called Einstein delay, as well as the deformation of space-time by the pulsar companions on the light travel path, the so-called Shapiro delay, and the aberration of the direction of the radio beam. All are calculated at 1PN order.

The delays due to interstellar medium propagation described by the DM as well as the Solar System counterparts of the previously mentioned delays are calculated by the commonly used software TEMPO2 (Edwards et al. 2006; Hobbs et al. 2006) which NUTIMO integrates as an external library. A thorough description of the timing model can be found in Chap. 5 of Voisin (2017).

3.1. Parametrisation of the problem

In total, the model must include at least 27 parameters (respectively 30 in the so-called secondary model when the 3 1PN strong-field parameters are included). One of them is a ToA uncertainty scale factor (called EFAC in the pulsar timing literature) which quantifies our ignorance of unmodelled systematic effects. The other 26 parameters (resp. 29) can be grouped into four categories:

• pulsar rotation: pulsar spin frequency and its derivative;

• orbital dynamics: six parameters for the inner-binary orbit, six parameters for the outer-WD orbit, three masses, one SEP violation parameter (resp. one SEP violation parameter and three 1PN strong-field parameters);

• astrometry: three position and three proper motion parameters;

• radio propagation: DM and DM derivative.

Each category is essentially uncorrelated with the others (see Fig. 2). The first two are specific to the triple-system problem and we shall discuss them in some details. On the other hand, the astrometric parameters, position and proper motion, and DM are treated using a standard approach and we refer the interested reader to Edwards et al. (2006), for example.

Fig. 2. Corner plot of the correlations between the fitted parameters of Table 1 and their marginalised distribution (diagonal histograms) sampled by our MCMC. A high-resolution version is available as online supplementary material.

The intrinsic pulsar parameters are its spin frequency f and spin-frequency derivative f′ taken at the reference epoch T_ref. These parameters need to be re-scaled to avoid non-linear correlations with astrometric parameters due to the Shklovskii delay (see Chap. 5 of Voisin 2017 for details about the delays). This is a common practice in pulsar timing. In addition, we re-scaled the spin frequency to include the linear effect (that is, proportional to time) of the Einstein delay, which is approximated by the second term of Eq. (5) below. In usual pulsar-timing models, the term of the Einstein delay responsible for a linear increase of the delay with time can be calculated exactly and removed from the timing model since its effect is strictly impossible to separate from a re-scaling of f. However, because here we calculate numerically the delay, we can only estimate the linear drift using the initial parameters. As a result of these re-scaling, the effective fit parameters $\overline{f}$ and ${\overline{f}}^{\prime }$ are connected to the intrinsic parameters f and f′ by

$$\begin{array}{cc}& \overline{f}=f[1-\frac{v^2+U_{\mathrm{i}}+U_{\mathrm{o}}}{2c^2}-{\mu }_{\perp }^2\mathrm{d}\mathrm{\Delta }T(1-\frac{3}{2}\mathrm{\Delta }T{\mu }_{\mathrm{d}})]/c,\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\overline{f}}^{\prime }=f^{\prime }-f\mathrm{d}{\mu }_{\perp }^2(1-3\mathrm{\Delta }T{\mu }_{\mathrm{d}})/c,\hfill \end{array}$$(6)

where

$$\begin{array}{cc}& v=2\pi (\frac{a_b}{P_{\mathrm{O}}}+\frac{a_{\mathrm{p}}}{P_{\mathrm{I}}}),\hfill \end{array}$$(7)

$$\begin{array}{cc}& U_{\mathrm{i}}=\frac{G{m}_{\mathrm{i}}}{a_{\mathrm{p}}(1+m_{\mathrm{p}}/m_{\mathrm{i}})},\hfill \end{array}$$(8)

$$\begin{array}{cc}& U_{\mathrm{o}}=\frac{G{m}_{\mathrm{o}}}{a_b(1+(m_{\mathrm{p}}+m_{\mathrm{i}})/m_{\mathrm{o}})},\hfill \end{array}$$(9)

$$\begin{array}{cc}& \mathrm{\Delta }T=T_{\mathrm{ref}}-T_{\mathrm{pos}},\hfill \end{array}$$(10)

where the symbols correspond to those defined in Table 1. The use of the two re-scaled parameters above instead of the intrinsic ones has proven to be very effective in speeding up convergence in our MCMC.

The orbital parameters for a triple system are at most 3 × 6 + 3 = 21 that is, three position coordinates and three velocity coordinates per body plus the three masses. However we consider the system in the frame of its centre of mass which results in applying two vector relations to the initial velocities and positions such that the centre of mass is at rest at coordinates (0, 0, 0). These relations suppress six degrees of freedom, and we end up with fifteen independent orbital parameters. Note that, eventually, the six degrees of freedom corresponding to the centre of mass appear as the six astrometric parameters.

The orbital parametrisation uses the fact that the system, being hierarchical, can approximately be described by an inner Keplerian binary containing the pulsar and the inner white dwarf (WD) itself forming an outer binary with the outer WD (see Fig. 3). Thus, the usual Keplerian orbital elements can be used to describe the osculating Keplerian orbits to the actual trajectory of the pulsar and of the inner binary. For eccentricity, we use the Laplace-Lagrange parametrisation relevant for small eccentricities (Lange et al. 2001) which replaces e, ω, t_p, respectively eccentricity, longitude of periastron and time of periastron passage, by e cos ω, e sin ω, and t_asc. It is important to note that we define the transformation t_asc = t_p − ω/2πP, with P an orbital period.

Fig. 3. Sketch of the orbits of PSR J0337+1715 (not to scale). Note that the orbits are nearly circular and that the ellipsoidal shape of the orbits on the left-hand sketch arises purely from projection. The osculating orbits of the system can be parametrised by the Keplerian orbital elements of the pulsar within the inner binary and of the inner binary within the outer binary (see text). In particular, ap is the semi-major axis of the pulsar orbit within the inner binary, ab is the semi-major axis of the inner binary within the outer binary and eO its eccentricity. Are also shown the longitude of periastron the outer binary as well as the inner and outer inclinations with respect to the plane of the sky, iI and iO respectively. Note that for simplicity the sketch neglects the difference between the directions of the inner and outer lines of ascending nodes, δΩ, as it is in practice very small.

Similarly, we fitted for a_b sin i_O, a_b cos i_O for the outer orbit, where a_b is the semi-major axis of the inner binary’s orbit around the centre of mass of the system and i_O its inclination relative to the plane of the sky. For the inner binary we find it better to fit for a_p sin i_I and δi = i_I − i_O instead of a_p sin i_I and a_p cos i_I as this cancels several non-linear correlations in the fit. This is helped by the fact the two orbits are very nearly coplanar. In the same way, we fitted for the longitude on the plane of the sky of the outer orbit, Ω_O, as well as for the offset between inner and outer orbits, δΩ = Ω_I − Ω_O. We note that only the latter was reported in Archibald et al. (2018) while the former was considered impossible to constrain with the available data. Interestingly, we were able to constrain Ω_O in this work, perhaps thanks to our inclusion of Kopeikin’s delays.

The inner binary mass m_b = m_i + m_p as well as the outer WD mass are derived using Kepler’s third law in the inner and outer binary respectively. The pulsar mass and the inner WD mass are derived from m_b using the mass ratio m_i/m_p which is also part of the fit. In addition, we use the post-Keplerian orbital elements of Damour & Deruelle (1985) which incorporate the 1PN corrections for relativistic binaries. Using Damour & Deruelle (1985), one maps the orbital elements to the position and velocity of the pulsar relative to the inner-binary centre-of-mass (r_p/b, v_p/b) and to the inner-binary centre-of-mass position and velocity (r_b, v_b). One can then find the position and velocity of the pulsar relative to the centre of mass of the system, r_p = r_p/b + r_b and v_i = v_i/b + v_b⁸. The position and velocity of the inner WD, (r_i, v_i), can then be derived by solving the equations of conservation of momentum and centre-of-mass position, (A.4) and (A.5) respectively, with P = (ℋ/c²)v_b and X = r_b. Finally, one solves P = 0 and X = 0 for the outer WD position and velocity, (r_o, v_o).

3.2. Model accuracy

The main signature of a SEP violation in our pulsar-timing experiment would be a residual signal at the frequency 2f_i − f_o (Archibald et al. 2018), where f_i, o are the inner and outer orbital frequency respectively (see also Fig. 4). According to linear theory, the effect of a violation of the SEP at Newtonian order primarily results in a sinusoidal variation of the separation of the inner binary with frequency f_i − f_o (Nordtvedt 1968). However, here we measure the distance projected along the line of sight between the observer and the pulsar and this distance is modulated by the orbit of the inner binary with the outer WD. It follows that the main net effect is a modulation at 2f_i − f_o as originally pointed out by Archibald et al. (2018). Using the Newtonian-order linear theory of Nordtvedt (1968), one can show that |Δ|∼10⁻⁶ translates into a signal amplitude of ∼0.1 s via the Rœmer delay. However, the magnitude that can effectively be detected in the fit residuals appears to be reduced to only ∼50 ns (Archibald et al. 2018) due to the effect of the many strong correlations of the Δ parameter with the other orbital parameters (see Fig. 2). Therefore, we aim in this work to achieve nanosecond accuracy within our model. This level of accuracy is compatible with the level aimed at by TEMPO2 (Edwards et al. 2006; Hobbs et al. 2006).

Fig. 4. Lomb-Scargle periodogram of the timing residuals for the solution of Table 1. The periodogram is sampled at fm/5 where fm ≃ 2268 days−1 is the inverse of the full time span. From left to right vertical lines show the frequencies of the red-noise component (fRN ≃ 1650 day−1, black), the Earth orbital period (fE, green), the outer-orbit orbital period (fO, orange), the inner-binary orbital period (fI, red) and the SEP violation signature (2fI − fO, purple).

There are essentially three types of inaccuracies that can affect the output of our model:numerical round-off errors, post-Newtonian truncation, interpolation precision. The first one, numerical round-off errors, is expected to grow with the number of floating-point operations performed to obtain the result. The main source of operations is the numerical integration of the equations of motion (A.2). The integration is performed using the Bulirsch-Stoer scheme (Stoer & Bulirsch 2011) implemented in the Odeint module (Ahnert et al. 2011) of the Boost library⁹. In addition our numerical model relies on 80 bit floating point numbers (long double in C) throughout. To assess the effect of numerical round-off errors, we use the model to generate fake times of arrival from parameters that fit the data from PSR J0337+1715. The fake times of arrival are the theoretical times of arrival that are closest to the actual measurements, and therefore only differ from those by a few microseconds at most. We feed this mock data set back into our model such that the residuals should be exactly zero in absence of numerical round-off errors. In practice we observe round-off errors at the level of 10⁻³ ns showing that numerical round-off errors are not an issue given our objective of a 1 ns accuracy. Note that this procedure does not assess any systematic inaccuracy due to the modelling or the numerical scheme themselves (see below) but it does account consistently for the entire chain of operations, including not only the numerical integration but also the Solar-system calculations done by TEMPO2 and the pulsar system delays. It is also conservative since the chain of operations is performed twice: once to create the fakes and once to analyse them.

The main source of systematic inaccuracy due to numerical approximations lies in the precision of the interpolations of the timing delays that are calculated in intermediate steps. We use a cubic-spline interpolation algorithm (Press 1996) for all our interpolations. There are two parameters than can be tuned: the number of interpolation points and the width of ‘margins’ at the beginning and the end of the interpolated range in order to avoid boundary effects. The latter need only be a few points in principle, however the former has a direct and opposite impact on accuracy and performance and therefore requires a trade-off. To determine the level of interpolation accuracy required we need to estimate what is the effect of a 1ns signal on the χ² value in order to make sure that this value is computed with the necessary accuracy. Let us assume that the difference between the data and the model is Δt_i + δ_i where the second term explicitly represent the contribution of a putative ∼1 ns signal, then the χ² can be expanded as follow,

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^N\frac{\mathrm{\Delta }{t}_i^2}{{\sigma }_i^2}+\sum _{i=1}^N\frac{\mathrm{\Delta }{t}_i{\delta }_i}{{\sigma }_i^2}+\sum _{i=1}^N\frac{{\delta }_i^2}{{\sigma }_i^2},}\end{array}$$(11)

where N is the number of data points. Now if N ≫ 1 and the number of fit parameters is ≪N, then for the best fit parameters Δt_i ∼ σ_i where {σ_i}_{i ∈ [1, N]} are the uncertainties which we take to be approximately equal to σ = 2 μs. In the present case N ∼ 10 000 for only 27 parameters. Thus, the first term in (11) approximates χ² ∼ N. The second term in Eq. (11) can be as large as Nδ/σ assuming that every term contributes positively. However it is also possible that the sum averages to zero if it alternates. The third term is of order ∼Nδ²/σ². Taking δ_i ∼ δ = 1 ns we see that the second and third term of equation (11) are respectively ≲5 × 10⁻⁴χ² and ∼2.5 × 10⁻⁷χ². We retain the last estimate as a conservative level of accuracy to achieve. To do so we increased exponentially the number of interpolation points until the relative variation of the χ² between two increments is smaller than 2.5 × 10⁻⁷.

The last source of inaccuracy, post-Newtonian truncation, is intrinsic to the theoretical framework used. Indeed, although the equations of motion and the various conserved quantities of Appendix A all consistently derive from the same Lagrangian and can therefore be exactly verified in principle, the method of calculus by successive approximation does not in practice achieve that result. Indeed, since the Lagrangian itself is an approximation to order v²/c² ∼ 5 × 10⁻⁷ of a more general theory, there is no physical justification for conserving in the subsequent derivation any term of higher order. It follows that the equations of motion and the corresponding conserved quantities are only accurate to first order (1PN) and that systematic ‘residuals’ of order v⁴/c⁴ ∼ 2 × 10⁻¹³ (2PN) are present in the equations themselves. As a result, we see on Fig. 5 that the energy of the system is conserved up to systematic oscillations at the orbital frequencies accompanied with a linear drift at a level consistent with the neglected 2PN terms, and numerical noise does not play any significant role. More interesting regarding the timing accuracy is to look at the conservation of the centre-of-mass position. Indeed, a shift in this position immediately transforms into a geometric delay. We find that, due to the fact that the neglected 2PN terms in the equations of motion do not necessarily generate residuals which average to zero, the two successive integrations $\ddot{x}\to \dot{x}\to x$ leading to the centre-of-mass motion create a quadratic drift that increases over time. Through the time span of our observations this results into a drift of less than 10 m, namely about 30 ns in terms of geometric delay. Such a quadratic drift can undoubtedly be entirely absorbed in the spin frequency and spin-frequency derivative as well as by the astrometric parameters when fitting the data. For example, the linear drift reported on Fig. 6 would bias the spin frequency by ∼10⁻¹⁴ Hz, much less than the nevertheless very tight uncertainty on this parameter. Therefore we conclude that the quadratic drift can only result in a negligible bias of a few parameters which is why we subtract this component with a linear-least-square fit on Fig. 6. The residuals show that the systematic oscillatory 2PN motion of the centre of mass does not exceed ∼0.1 m, or 0.3 ns in terms of geometric delay, which is well within our tolerance.

Fig. 5. Relative error in energy conservation during the numerical integration of the equations of motion over the 2268 day span of our data. The envelope of the signal oscillates at the outer binary frequency while a zoom would show a fast oscillation at the inner binary orbital frequency.

Fig. 6. Three components of the centre-of-mass position variation during the integration of the equations of motion over the 2268 day span of our data. A quadratic component has been fitted out (see text and legend) which leaves only the oscillatory components. The main visible oscillation is at the frequency of the outer binary while a zoom would show a fast oscillation at the inner binary orbital frequency.

In practice, the largest systematic errors may come from unmodelled effects. In particular, gravitational wave damping in the inner binary should account for a few nanoseconds after 5 years. Another effect that might become important for high-precision timing over time is the effect of the gravitational quadrupole moment of the inner white dwarf. Indeed this star should be slightly deformed by the tidal field of the neutron star and by its spin which would lead to a slight correction to the orbital precession rate.

4. Bayesian analysis results

Our main goal in this work has been to get a Bayesian estimate of the uncertainties on each of the parameters of the problem, or in other words to estimate the posterior probability density function (PDF) of the parameters θ belonging to model M given our data D using Bayes’ rule,

$$\begin{array}{c}\hfill P(\mathit{\theta },M|D)=\frac{P(D|\mathit{\theta },M)P(\mathit{\theta },M)}{P(D)}·\end{array}$$(12)

The prior function, P(θ, M), was chosen flat except for astrometric parameters that benefited from prior knowledge of position and angular proper motion from the Gaia mission DR2 (Lindegren et al. 2018), of distance from photometric observations of the inner white dwarf (Ransom et al. 2014) and radial velocity from optical spectroscopy of the same star (Kaplan et al. 2014). The Gaia DR2 catalogue does not model orbital motion which may then contaminate both position and proper motion. In the present case, given the distance of the source the magnitude of orbital motion is similar to the uncertainties reported by Gaia DR2. In order to account for potential systematic errors we have multiplied by two these uncertainties before using them as standard deviations of our Gaussian priors (see Table 2). We have also applied a factor of two to the uncertainty on the photometric distance reported in Ransom et al. (2014) in order to account for potential systematic effects that would bias this measurement. For instance an inaccurate spectroscopic estimate of surface gravity (the ‘high log g problem’ in low-mass white dwarfs, see Tremblay et al. 2015 and references therein) would in turn bias the stellar radius estimate and therefore the absolute magnitude of the star. It is worth pointing out that the two commonly used free-electron density models for the Galaxy, NE2001 (Cordes & Lazio 2002) and YMW16 (Yao et al. 2017), both predict a distance of about 800 pc which is significantly smaller that reported in Tables 1–2. This indicates that the electron density for the given Galactic height (z = −690(40) pc) is overestimated. All priors are summarised in Table 2. Let us note that our fit for the radial proper motion μ_d is unconstraining as the uncertainties reported in Table 1 match the radial velocity prior of Table 2. The uncertainties of all the other fitted quantities are improved with respect to their prior.

The high dimensionality of the PDF together with the necessity to integrate numerically the equations of motion makes the problem computationally challenging. However our C++ code is able to calculate one PDF value in less than 10 s on a last-generation laptop, which made it possible to sample the PDF on a medium-size computer cluster. The sampling was achieved using a home-made implementation of the affine-invariant Markov-chain Monte Carlo (MCMC) of Goodman & Weare (2010) parallelised with the scheme of Foreman-Mackey et al. (2013). The advantage of this algorithm is to be efficient in high dimensionality (Allison & Dunkley 2014) and insensitive to any level of linear correlations between the parameters. This is particularly important as we found Δ to be highly correlated with many orbital parameters (see Fig. 2). However, we also found that non-linear ‘correlations’ between parameters were preventing convergence within a reasonable time, which was solved by appropriate re-parametrisation (see Sect. 3.1). Convergence was evaluated by requiring that fluctuations of the mean and standard-deviation estimators be smaller that 6% of the full-chain standard deviation for each parameter (see e.g. Dunkley et al. 2005, Sect. 4.1 and chain plot in supplementary online material). We noticed that standard deviations sometimes converge later than means, particularly for Δ, confirming the importance of monitoring both indicators to ensure reliable uncertainties.

Due to its very low ecliptic latitude, ∼2 deg, the timing of PSR J0337+1715 is potentially sensitive to a range of effects occurring in the Solar-system. In particular, we detected in preliminary runs a slight increase in timing residuals of the order of 1 μs when the pulsar was within 3 deg of the Sun. We attributed this increase to the inaccuracy of the Solar-wind electron density model used by TEMPO2 to calculate the related DM. We mitigated this effect by removing all the ToAs taken within 5 deg of the Sun. Moreover, our periodogram shows a secondary ∼0.2 μs component close to the Earth orbital frequency, sign of possible extra inaccuracies in the Solar-wind model or in the Solar-system ephemerides. This is likely to affect outer-orbit parameters since this period is close to 1 year but such a correlation can only widen posterior uncertainties.

Two models M were tested. Our main model includes only Δ as a free parameter while our secondary model includes the three additional 1PN-strong-field parameters, ${\overline{\beta }}_{\mathrm{p}}$, ${\overline{\beta }}_{\mathrm{pp}}$, ${\overline{\beta }}^{\mathrm{p}}$, yielding the following 95% C.L. constraints for them:

$$\begin{array}{cc}& \mathrm{\Delta }=1.1_{-4.6}^{+4.8}\times {10}^{-6},\hfill \end{array}$$(13)

$$\begin{array}{cc}& {\overline{\beta }}_{\mathrm{p}}=-2_{-14}^{+13},\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\overline{\beta }}_{\mathrm{pp}}=-0.1_{-1.0}^{+0.9},\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\overline{\beta }}^{\mathrm{p}}=+0.8_{-6.6}^{+7.2}.\hfill \end{array}$$(16)

In regard to binary-pulsar tests (see Sect. 5), the above results on the three $\overline{\beta }$ parameters are unconstraining. We used this prior knowledge to run our main model with ${\overline{\beta }}_{\mathrm{p}}={\overline{\beta }}^{\mathrm{p}}={\overline{\beta }}_{\mathrm{pp}}=0$ and obtain our primary SEP limit

$$\begin{array}{c}\hfill \mathrm{\Delta }=(+0.5\pm 1.8)\times {10}^{-6}(95\%\mathrm{C}.\mathrm{L}.),\end{array}$$(17)

which translates into |Δ| < 2.05 × 10⁻⁶ at 95% C.L. (see Fig. 7). The full result of the main model is reported in Table 1. Note that ∼8% of the reported uncertainties are due to unaccounted systematics absorbed in the EFAC parameter (see also Sect. 4.2). The wider uncertainty obtained in the secondary run is due to large correlations with the three additional parameters.

Fig. 7. Left-hand side: marginalised posterior probability distribution of the SEP violation parameter Δ sampled by MCMC and normal law with the same mean and standard deviation. The upper axis gives the mean value and the boundaries of the 95% confidence region. Right-hand side: distribution of distance to GR derived from the left-hand-side distribution.

4.1. MCMC run and convergence

The affine-invariant algorithm of Goodman & Weare (2010) requires to move N walkers together at each iteration. The gist of this algorithm is that the walkers within the set are not independent from each other while the set as a whole constitutes a single effective walker in the Markov process sense, namely that it depends only on its previous state. Individual moves within the set are informed by the positions of other walkers in a way that renders the algorithm rigorously immune to any linear correlation, or any affine parameter transformation. However it might be sensitive to correlation of a higher degree, or to non-convexity of the posterior isosurfaces. Therefore, with this algorithm one should take care of removing as much as possible any non-linear correlations by choosing an appropriate parameter set (see Sect. 3.1) but very large linear correlations, as can be seen in Fig. 2, are well resolved by the algorithm.

The authors of Goodman & Weare (2010), Foreman-Mackey et al. (2013) recommend to choose a number of walkers within the set much larger than the number of parameters. In the present case we chose to use 288 walkers per set. The only other tunnable parameter is the unique parameter of the proposal function, a, which controls the size of the steps that can be attempted. The authors of Goodman & Weare (2010), Foreman-Mackey et al. (2013) suggest the value a = 2 in order to keep an acceptance fraction ∼0.4. We found that, as the chain was getting closer to convergence the acceptance fraction could drop dramatically, sign of non-linear correlations or non-convexity. This drop was largely mitigated by adopting the final parameter set of Sect. 3.1, but we still had to choose a smaller step-size parameter a = 1.6 to keep the acceptance fraction close to the level mentioned above.

Due to the large parameter space and the computing time needed to compute one χ² (about 10 s) we parallelised the MCMC code using the scheme of Foreman-Mackey et al. (2013). This allowed us to run the MCMC using 144 cores of the meso-scale MesoPSL cluster (see acknowledgements) each calculating 2 walkers (which is the minimum number of walkers per core given the algorithm used). A few 10 000 steps were typically necessary to reach convergence. However, it is important to quantitatively estimate convergence as one cannot afford to run the MCMC for an unnecessarily large number of iterations. Beyond visual inspection of the parameter chains which allows to discard any obvious burn-in phase, we also monitored the autocorrelation time of each chain (Goodman & Weare 2010). However, we found that the most constraining convergence estimator was to monitor the variation of the mean and standard deviation of each parameter chain. Formally, one needs to compute the value of the following estimator on each chain (Dunkley et al. 2005),

$$\begin{array}{c}\hfill {\widehat{r}}_{\widehat{o}=\widehat{m},\widehat{\sigma }}=\frac{\widehat{\sigma }\left(\widehat{o}\right)}{\sigma },\end{array}$$(18)

where $\widehat{\sigma }$ is a standard-deviation estimator, $\widehat{o}$ is a statistical estimator which here is either the mean $\widehat{m}$ or the standard deviation $\widehat{\sigma }$, and σ is the standard deviation estimated on the entire length of the chain. In practice, we recorded the 288 walkers every 5 iterations and use the last 69 408 recorded elements (241 independent sets of 288 walkers). The standard deviation of $\widehat{o}$ was estimated by applying $\widehat{o}$ on 20 sub-samples of the chain and then estimating the standard deviation of the set of the values obtained. The chain was considered converged if both ${\widehat{r}}_{\widehat{m}}$ and ${\widehat{r}}_{\widehat{\sigma }}$ estimated values are smaller than 0.06. A situation we observed is when the latter keeps varying significantly while the former is stable and under the cutoff. In other words, the chain widens with constant mean, rendering a mean-based convergence estimator uninformative and possibly leading to an underestimation of the parameter uncertainties.

4.2. Analysis of the residuals

We have assessed the robustness of our fit by evaluating the presence of systematic components in the residuals (Fig. 8). As it appears from Fig. 9, no significant structure is present at either the inner or outer orbital period, nor at the Earth orbital period notwithstanding the sharp cut around the phase of closest approach to the Sun made to prevent any bias induced by unmodelled DM contributions.

Fig. 8. Residuals of the times of arrival using the mean parameters reported in Table 1.

Fig. 9. Residuals of the times of arrival corresponding to the parameters of Table 1 plotted versus the inner-binary orbital phase (top), the outer-binary orbital phase (middle) and the Earth orbital phase (bottom). The red vertical bars of the first two plots show the phase where the pulsar (resp. the inner binary) is closer to the Solar System. The red vertical bar on the bottom plot, shows the Earth orbital phase where the line of sight to the pulsar passes closest to the Sun. We have removed every ToA when the pulsar is less than 5 deg from the Sun, which explains the gap around that particular phase.

In order to estimate more thoroughly the presence of systematic modulations we produced a Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) of the same residuals, Fig. 4. It appears that indeed there is no significant power at the frequencies mentioned above, except for a tentative signal ≲0.2 μs at the Earth orbital frequency. Because of the proximity of this frequency with the outer-binary orbital frequency this might lead to correlate Solar-system and outer-binary effects and therefore enlarge uncertainties related to the parameters involved.

The dominant component is the low-frequency peak and its subsequent harmonics which we interpret as time-correlated red noise. A number of causes have been proposed in the literature. The main ones are the intrinsic spin frequency noise (Shannon & Cordes 2010; Melatos & Link 2014) or magnetospheric fluctuations (Lyne et al. 2010). It has also been proposed that asteroid belts could result in apparent timing noise (Shannon et al. 2013). Propagation effects under the form of long-term variations of the dispersion measure along the line of sight due to turbulence in the interstellar medium could be the cause in some cases (Keith et al. 2013). We have tried to split the time span of our observations into several intervals with different DM values, but the fit was consistent with an absence of variation discarding this explanation. Red noise can also have a local cause, such as irregularities of the terrestrial time realisation used to time the pulsar (Hobbs et al. 2020, 2012) or inaccuracies in the planetary masses used to generate the Solar System ephemeris (Champion et al. 2010; Caballero et al. 2018). However, we use the 2016 realisation of the BIPM terrestrial time which Hobbs et al. (2012) has confirmed as suitable for precision pulsar timing. Besides, if the red noise was caused by any inaccuracy in planetary masses, then the signature would be at the orbital frequency of the responsible planet (Champion et al. 2010). The only orbital period in the Solar System that approaches the 1650 days of the red-noise signal is the orbital period of the dwarf planet Ceres. However the uncertainty on its mass in the NASA JPL ephemeris DE430 (Folkner et al. 2014) we use in this work is too small to explain a signal of that magnitude.

Thus, there is no deterministic model that can be fitted to that component, but since its frequency is much lower than any other in the system it is unlikely to bias the parameters. However, it results in increasing the scale of the ToA uncertainties (EFAC parameter in Table 1) in order to accommodate this systematic effect into a reduced χ² equal to unity. Were the PDF perfectly Gaussian, that would result into posterior uncertainties increased in exactly the same proportion, which we can here estimate at ∼8%. Therefore, our posterior uncertainties should be seen as upper limits. Interestingly, because the analysis of Archibald et al. (2018) focuses on a specific frequency signature for the SEP, the result quoted in that work should be seen as a lower limit on the actual uncertainty in the sense of the above discussion.

Finally, the prominent peak at 1 day⁻¹ and its harmonics on the periodogram of Fig. 4 result from the convolution of the red noise component with the observing window functions of our data. Indeed, due to its meridian configuration, the Nançay radio telescope can only observe the same object during ∼1 h windows separated by an integer number of sidereal day. Due to the proximity of the inner orbital frequency with the observing frequency of ∼1 day⁻¹ one might be concerned with a potential bias. However we have checked that the Fourier transform of a comb of 1 h window functions results in sharp narrow peaks whose width does not exceed a few percents of the daily frequency and therefore cannot significantly bias signals at the orbital frequencies.

5. Theory dependent tests

The parametrised post-Newtonian (PPN) formalism (see e.g. Will 2018a), with its ten theory-independent parameters, has proven to be a powerful tool to quantify and compare tests of GR and its alternatives in the weak-field environment of the Solar System. Unfortunately, there is no such framework that generically extends beyond the weak field approximation of the PPN formalism and therefore is able to cover gravity experiments with strongly self-gravitating bodies, like the one in this paper. One reason is that, unlike in the Solar System, the treatment of the motion of a neutron star in an external gravitational field requires the full complexity of a gravity theory (Will 2018b). To put the UFF test conducted in this paper into context with other experiments, in particular Solar System tests and gravitational wave tests with binary pulsars, it has been suggested to use theory-dependent frameworks (see e.g. Damour 2009). Scalar-tensor theories of gravity have turned out to be particularly useful for this purpose. Apart from being well motivated and well studied alternatives to GR (Fujii & Maeda 2007), they show a rich phenomenology in their deviations from GR, including prominent effects related to the non-linear strong-field regime of neutron stars (see e.g. Damour & Esposito-Farèse 1993, 1996).

In this paper, as a theory-dependent framework we use the class of Bergmann-Wagoner theories. Bergmann-Wagoner theories represent the most general scalar–tensor theories with one scalar field that are at most quadratic in the derivatives of the fields (Will 2018a). Quite a number of well known scalar-tensor theories belong to this class, like Jordan-Fierz-Brans-Dicke (JFBD) gravity (Jordan 1955; Fierz 1956; Brans & Dicke 1961), DEF gravity (Damour & Esposito-Farèse 1993), MO gravity (Mendes & Ortiz 2016), f(R) gravity (De Felice & Tsujikawa 2010), and massive Brans-Dicke gravity (Alsing et al. 2012). Bergmann-Wagoner theories form a subclass of the class of Horndeski theories (Horndeski 1974), which is the most general class of mono-scalar-tensor theories in four dimensions whose Lagrangian leads to second order field equations.

In the following we interpret our limits of Sect. 4 in two different approaches within the class of Bergmann-Wagoner theories. In the first approach we remain (mostly) generic, in the sense that make as few assumptions as possible concerning the details of the theory. In the second approach, we pick a specific two-parameter scalar-tensor theory, that is, DEF gravity. For this two-parameter class of theories we can then explicitly calculate the properties of neutron stars for different equations of state (EOS) and derive constraints on the two-dimensional theory space from the limits presented here.

5.1. Generic tests within Bergmann-Wagoner scalar-tensor gravity

In Bergmann-Wagoner theories, the field equations for the (physical) metric g_μν and the scalar field ϕ are a result of the following action

$$\begin{array}{cc}\hfill S=& \frac{1}{16\pi G_0}{\displaystyle \int \sqrt{-g}{d}^{4}x(\varphi R-\frac{\omega (\varphi )}{\varphi }{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi -U(\varphi ))}\hfill \\ \hfill & +S_{\mathrm{mat}}[{\psi }_{\mathrm{mat}}^{\mathrm{A}},g_{\mu \nu }],\hfill \end{array}$$(19)

where G₀ is the fundamental (‘bare’) gravitational constant, g the determinant of g_μν, R the curvature scalar, ω(ϕ) is the coupling function, and U(ϕ) the scalar potential. The physical (Newtonian) gravitational constant, as measured in a Cavendish-type experiment, is given by

$$\begin{array}{c}\hfill G_{\mathrm{N}}=\frac{G_0}{{\varphi }_0(1-\zeta )},\end{array}$$(20)

where ϕ₀ denotes the cosmological background field and ζ ≡ 1/(2ω(ϕ₀)+4). S_mat is the action of the matter fields ${\psi }_{\mathrm{mat}}^{\mathrm{A}}$, which couple universally to the spacetime metric g_μν. For our discussion, we assume that U(ϕ) can be neglected on the scale of the triple system, that is, $U^{″}(\varphi )\ll 1/a_b^2$. In terms of a massive scalar field, this means that we assume that the Compton wavelength is much larger than the extension of the system (see Seymour & Yagi 2019, on how J0337+1715 can be used to constrain massive scalar fields).

The effective gravitational constant that enters the N-body Lagrangian is given by

$$\begin{array}{c}\hfill G_{\mathit{ab}}=G_{\mathrm{N}}[1-2\zeta (s_a+s_b-2s_a{s}_b)],\end{array}$$(21)

where the sensitivity

$$\begin{array}{c}\hfill s_a\equiv {\left(\frac{\mathrm{d}ln{m}_a(\varphi )}{\mathrm{d}ln\varphi }\right)}_{{\varphi }_0}\end{array}$$(22)

accounts for the dependence of each body on a change in the ambient scalar field, while the number of baryons remains fixed (Will 2018a). For neutron stars, s_a depends on the EOS. It is typically of the order of 0.1 but, depending on the details of ω(ϕ), its (absolute) value can be very large, as we discuss further below. For the relative GWEP parameter one finds

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathit{ab}}=-2\zeta (s_a+s_b-2s_a{s}_b).\end{array}$$(23)

Because of the product s_as_b, in general it is not possible to interpret the quasi-Newtonian equations of motion in terms of inertial and gravitational masses of the individual bodies in an N-body system (Will 2018a). For weakly self-gravitating bodies the sensitivity s_a is simply related to the fractional gravitational binding energy ε_grav, a via

$$\begin{array}{c}\hfill s_a=-(1+2\lambda ){\epsilon }_{\mathrm{grav},a}+\mathcal{O}\left({\epsilon }_{\mathrm{grav},a}^2\right),\end{array}$$(24)

where λ ≡ ϕ₀ω′(ϕ₀)ζ²(1 − ζ)⁻¹. The two Eddington parameters of the PPN formalism are given by

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{PPN}}=1-2\zeta ,{\beta }_{\mathrm{PPN}}=1+\zeta \lambda ,\end{array}$$(25)

and the Nordtvedt parameter of Eq. (4) reads

$$\begin{array}{c}\hfill \eta =4{\beta }_{\mathrm{PPN}}-{\gamma }_{\mathrm{PPN}}-3=2\zeta (1+2\lambda ).\end{array}$$(26)

For the inner and outer white dwarf we have ε_grav, i ≃ −1.8 × 10⁻⁵ and ε_grav, o ≃ −4.8 × 10⁻⁵, respectively. Consequently G_io ≃ G_N for the interaction between the two white dwarfs, and G_pi ≃ G_po ≃ G_N(1 − 2ζs_p) for the interaction between the pulsar and the white dwarfs. Hence, our result for Δ in Table 1 leads to a direct constraint for

$$\begin{array}{c}\hfill \zeta s_{\mathrm{p}}\simeq -\frac{\mathrm{\Delta }}{2}=(-0.2\pm 0.9)\times {10}^{-6}(95\%\mathrm{C}.\mathrm{L}.),\end{array}$$(27)

where Δ ≡ Δ_pb(s_b = 0). The above limit can be considered as generic within the family of Bergmann-Wagoner theories of gravity, in the sense that it does not require a specification of the coupling function ω(ϕ). Later, we use this limit to impose constraints on the parameter space of a specific two parameter family of Bergmann-Wagoner theories. Before that, we need to discuss the strong-field modifications at the post-Newtonian level of the 3-body dynamics.

5.1.1. First post-Newtonian contributions

At the first post-Newtonian (1PN) level (order v²/c²), the PPN parameters ${\overline{\beta }}_{\mathrm{PPN}}\equiv {\beta }_{\mathrm{PPN}}-1$ and ${\overline{\gamma }}_{\mathrm{PPN}}\equiv {\gamma }_{\mathrm{PPN}}-1$ need to be replaced by the body-dependent quantities ${\overline{\gamma }}_{\mathit{ab}}$ and ${\overline{\beta }}_{\mathit{bc}}^a$ (see Appendix A for details). These strong-field generalisations of the PPN Eddington parameters depend on the sensitivities and their derivatives of the bodies in a system. For the detailed expressions, we refer the reader to Will (1993) and Damour & Esposito-Farèse (1992). The latter uses the so-called Einstein-frame representation and gives these terms for multi-scalar-tensor theories. More generally, in the triple system of PSR J0337+1715, where two of the bodies are weakly self gravitating, one finds for the twelve 1PN-strong-field parameters, to good approximation,

$$\begin{array}{cc}& {\overline{\gamma }}_{\mathrm{io}}\simeq {\overline{\gamma }}_{\mathrm{PPN}},\hfill \end{array}$$(28)

$$\begin{array}{cc}& {\overline{\gamma }}_{\mathrm{pi}}\simeq {\overline{\gamma }}_{\mathrm{po}}\simeq {\overline{\gamma }}_{\mathrm{p}b}(s_b=0)\equiv {\overline{\gamma }}_{\mathrm{p}},\hfill \end{array}$$(29)

$$\begin{array}{cc}& {\overline{\beta }}_{\mathrm{oo}}^{\mathrm{i}}\simeq {\overline{\beta }}_{\mathrm{ii}}^{\mathrm{o}}\simeq {\overline{\beta }}_{\mathrm{PPN}},\hfill \end{array}$$(30)

$$\begin{array}{cc}& {\overline{\beta }}_{\mathrm{po}}^{\mathrm{i}}\simeq {\overline{\beta }}_{\mathrm{pi}}^{\mathrm{o}}\simeq {\overline{\beta }}_{\mathrm{p}c}^a(s_a=s_c=0)\equiv {\overline{\beta }}_{\mathrm{p}},\hfill \end{array}$$(31)

$$\begin{array}{cc}& {\overline{\beta }}_{\mathrm{pp}}^{\mathrm{i}}\simeq {\overline{\beta }}_{\mathrm{pp}}^{\mathrm{o}}\simeq {\overline{\beta }}_{\mathrm{pp}}^a(s_a=0)\equiv {\overline{\beta }}_{\mathrm{pp}},\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\overline{\beta }}_{\mathrm{ii}}^{\mathrm{p}}\simeq {\overline{\beta }}_{\mathrm{io}}^{\mathrm{p}}\simeq {\overline{\beta }}_{\mathrm{oo}}^{\mathrm{p}}\simeq {\overline{\beta }}_{\mathit{bc}}^{\mathrm{p}}(s_b=s_c=0)\equiv {\overline{\beta }}^{\mathrm{p}},\hfill \end{array}$$(33)

leaving us with six different parameters at the 1PN level of the modified Einstein-Infeld-Hoffmann equations of motion, instead of the two in the weak field limit. Note the following symmetries: ${\overline{\gamma }}_{\mathit{ab}}={\overline{\gamma }}_{\mathit{ba}}$ and ${\overline{\beta }}_{\mathit{bc}}^a={\overline{\beta }}_{\mathit{cb}}^a$.

At this stage, we can further reduce the number of 1PN parameters, without making more detailed assumptions about the theory, for instance about ω(ϕ). The tight limits on ${\overline{\beta }}_{\mathrm{PPN}}$ and ${\overline{\gamma }}_{\mathrm{PPN}}$ from Solar System tests ∼10⁻⁵ (Will 2018a), directly put tight constraints on two of the six 1PN parameters. Furthermore, ${\overline{\gamma }}_{\mathrm{p}}$ only depends on terms proportional to ζ and ζs_p, the first being constrained to ∼10⁻⁵ by Cassini (Bertotti et al. 2003) and the latter to ∼10⁻⁶ already by the Newtonian-level dynamics of the triple system (cf. Eq. (27); see also limit (13)). Consequently, without loss of generality, ${\overline{\gamma }}_{\mathit{ab}}$, ${\overline{\beta }}_{\mathrm{oo}}^{\mathrm{i}}$, and ${\overline{\beta }}_{\mathrm{ii}}^{\mathrm{o}}$ can be ignored in a self-consistent gravity test with the PSR J0337+1715 system. Besides the Δ at the Newtonian level, we are left with the 1PN strong-field parameters ${\overline{\beta }}_{\mathrm{p}}$, ${\overline{\beta }}_{\mathrm{pp}}$, and ${\overline{\beta }}^{\mathrm{p}}$. These three parameters cannot be constrained without further assumptions, as we discuss below. Hence we have implemented a model based on deviations from GR parametrised by $(\mathrm{\Delta },{\overline{\beta }}_{\mathrm{p}},{\overline{\beta }}_{\mathrm{pp}},{\overline{\beta }}^{\mathrm{p}})$, which is called secondary model in Sects. 3 and 4. Our analysis based on this model leads to the generic limits (13)–(16).

In our generic approach, the parameters ${\overline{\beta }}_{\mathrm{p}}$, ${\overline{\beta }}_{\mathrm{pp}}$ can a-priori only be constrained if we make further assumptions and apply existing constraints from binary-pulsar systems. The reason is as follows. The three $\overline{\beta }$ parameters have terms, which are proportional to λζ, ${\zeta }^2{s}_{\text{p}}^2$, λζs_p, $\lambda \zeta s_{\text{p}}^2$, and ${\zeta }^2{s}_{\text{p}}^{\prime }$, where

$$\begin{array}{c}\hfill s_a^{\prime }\equiv {\left(\frac{{\mathrm{d}}^{2}ln{m}_a(\varphi )}{\mathrm{d}{(ln\varphi )}^2}\right)}_{{\varphi }_0}\end{array}$$(34)

(see e.g. Will 2018b, for details). Solar System constraints on β_PPN and η put tight constraints (∼10⁻⁵) on λζ, and ${\zeta }^2{s}_{\text{p}}^2$ is constrained to ∼10⁻¹² because of Eq. (27). However, the quantities λ, s_p, and $s_{\text{p}}^{\prime }$ are unconstrained by above considerations. Consequently, λζs_p, $\lambda \zeta s_{\text{p}}^2$, and ${\zeta }^2{s}_{\text{p}}^{\prime }$ are a-priori unconstrained. Below, in a more theory specific context, we discuss a situation in DEF gravity where $\zeta s_{\text{p}}^2$ can be of order unity, even if ζs_p ≪ 1 (spontaneous scalarisation of neutron stars). In such a highly non-linear strong field regime, s_p and $s_{\text{p}}^{\prime }$ can assume very large values. As a consequence, the $\overline{\beta }$ terms can become much larger than β_PPN, even under condition (27). To give an example, in the regime of spontaneous scalarisation in DEF gravity, one finds ${\overline{\beta }}_{\mathrm{pp}}/{\overline{\beta }}_{\mathrm{PPN}}\sim {\zeta }^{-1}$ while ${\overline{\beta }}_{\mathrm{pp}}$ remains practically unaffected when ζ → 0 (cf. Damour & Esposito-Farèse 1996).

If we restrict λ to values which are not very large, it can be shown that ${\overline{\beta }}_{\mathrm{p}}$ can be considered as small as well, since it only contains λζs_p as an a priori unconstrained term, where ζs_p has to be small according to Eq. (27).

Further constraints can come from binary pulsar experiments. In particular, dipolar-radiation tests in binary pulsars with spectroscopic white dwarfs (see e.g. Lazaridis et al. 2009; Freire et al. 2012b), can in principle provide generic constraints on $\zeta s_{\text{p}}^2$, although some additional assumptions are needed, for instance for U(ϕ) (cf. Alsing et al. 2012). If U(ϕ) = 0, the change in the orbital period of a pulsar-white dwarf binary, P_b, due to dipolar radiation damping is given by

$$\begin{array}{c}\hfill {\dot{P}}_b^{\mathrm{dipole}}\simeq -\frac{G_{\mathrm{N}}}{c^3}\frac{16{\pi }^2}{P_b}\frac{m_{\mathrm{p}}{m}_c}{m_{\mathrm{p}}+m_c}\frac{1+e^2/2}{{(1-e^2)}^{5/2}}\zeta s_{\mathrm{p}}^2\end{array}$$(35)

(see e.g. Eq. (12.32) in Will 2018a), where e denotes the orbital eccentricity, and m_p and m_c are the masses of pulsar and white-dwarf companion respectively. To apply constraints from other pulsar to the PSR J0337+1715 system it also requires some on the dependence of the sensitivity of the pulsar, s_p, on the pulsar mass. In the strong-field regime of neutron stars this dependence can be highly non-linear (Damour & Esposito-Farèse 1996; Shao et al. 2017). PSR J1738+0333 is a pulsar with a mass similar to PSR J0337+1715 (m_p ≃ 1.46 M_⊙). The dipolar radiation test by Freire et al. (2012b); Zhu et al. (2019) leads to

$$\begin{array}{c}\hfill \zeta s_{\mathrm{p}}^2=(-0.5\pm 1.7)\times {10}^{-6}(95\%\mathrm{C}.\mathrm{L}.).\end{array}$$(36)

As a result, ${\overline{\beta }}_{\mathrm{pp}}$ can, in general, also be assumed to be small in the PSR J0337+1715 system.

Imposing a generic constraint on ${\zeta }^2{s}_{\text{p}}^{\prime }$, and therefore on ${\overline{\beta }}^{\mathrm{p}}$, is somewhat less direct. $s_{\text{p}}^{\prime }$ enters, for instance, the precession of periastron, $\dot{\omega }$, which is particularly well tested – in combination with other post-Keplerian parameters – in eccentric short-orbital-period binary pulsars (Wex 2014; Will 2018a). However, only the so called Double Pulsar allows for a generic constraint on deviations from GR in $\dot{\omega }$, which is of the order of 10⁻³ (Kramer & Wex 2009). However, the masses of the Double Pulsar are significantly lower than the mass of PSR J0337+1715. Nevertheless, the general agreement of all these systems with GR at least suggests that ${\zeta }^2{s}_{\text{p}}^{\prime }$, and therefore ${\overline{\beta }}^{\mathrm{p}}$ can generally be assumed to be small as well. Moreover, given that the Cassini experiment already imposes ζ² ≲ 10⁻¹⁰, $s_{\text{p}}^{\prime }$ would have to assume quite extreme values to lead to a significant ${\overline{\beta }}^{\mathrm{p}}$, certainly in view of the (still) quite weak limit (16).

To summarise, under additional assumptions, which we consider as reasonable for most situations, all three strong-field $\overline{\beta }$ parameters are tightly constrained by a combination of Eq. (27) with constraints from binary pulsars experiments. The limits in Eqs. (14)–(16) are therefore generally not of particular interest, at least for constraining the class of scalar-tensor theories considered in this section.

5.2. EOS-agnostic constraints on Damour–Esposito-Farèse (DEF) gravity

In order to explicitly calculate s_p and $s_{\text{p}}^{\prime }$ and therefore Δ and the 1PN strong-field parameters that enter our equations of motion, one has to pick a specific theory of gravity, which we do in this subsection. In the quadratic mono-scalar tensor theory T₁(α₀, β₀) of Damour & Esposito-Farèse (1993), the coupling function in the Einstein frame is quadratic in the scalar field, meaning that the coupling strength between the scalar field and the trace of the stress-energy tensor becomes field dependent in a linear way. In the Jordan-frame representation with the physical metric g_μν, which we are using here, the coupling function then reads

$$\begin{array}{c}\hfill \omega (\varphi )=\frac{1}{2}(\frac{1}{{\alpha }_0^2-{\beta }_{0}ln\varphi }-3),\end{array}$$(37)

where ϕ₀ ≡ 1, without loss of generality (Will 2018a). Furthermore, one finds

$$\begin{array}{c}\hfill \zeta =\frac{{\alpha }_0^2}{1+{\alpha }_0^2},\lambda =\frac{{\beta }_0}{2(1+{\alpha }_0^2)}·\end{array}$$(38)

The tight constraints on ζ from the Cassini mission imply that ${\alpha }_0^2\lesssim {10}^{-5}$. Furthermore, from Eq. (26) one then finds for the Nordtvedt parameter

$$\begin{array}{c}\hfill \eta \simeq 2{\alpha }_0^2(1+{\beta }_0).\end{array}$$(39)

The sensitivity s_a of a weakly self-gravitating body, like the WD companions to J0337+1715, can be calculated according to Eq. (24):

$$\begin{array}{c}\hfill s_a\simeq -(1+{\beta }_0){\epsilon }_{\mathrm{grav},a}.\end{array}$$(40)

For neutron stars, the absolute value of the sensitivity can become very large if β₀ ≲ −4.5, a fact first discovered within T₁(α₀, β₀) gravity theories by Damour & Esposito-Farèse (1993), and generally referred to as ‘spontaneous scalarisation’. As a result, even for arbitrarily small α₀, the quantity α₀s_a ≃ ζ^1/2s_a remains at order unity¹⁰.

A special case of T₁(α₀, β₀) is JFBD gravity, for which β₀ = 0. In that case λ = 0, and the coupling function is a constant:

$$\begin{array}{c}\hfill \omega (\varphi )=\frac{1}{2{\alpha }_0^2}-\frac{3}{2}=\frac{1}{2\zeta }-2\equiv {\omega }_{\mathrm{BD}}>-\frac{3}{2},\end{array}$$(41)

which is called the Brans-Dicke parameter. For ω_BD → ∞, that is, α₀ and ζ → 0, JFBD gravity approaches GR. We obtain the most conservative limits on JFBD when using the stiffest EOS from our set of viable EOSs (see Fig. 10), that is, BSk22. For this EOS, in JFBD gravity, the sensitivity of PSR J0337+1715 has the value s_p = 0.149. Most importantly, for ζ ≲ 10², this value is practically independent of ζ (Shibata et al. 2014; Shao et al. 2017). Hence, Eq. (27) can directly be converted into limits on the coupling parameter:

$$\begin{array}{c}\hfill \zeta =(-1\pm 6)\times {10}^{-6}(95\%\mathrm{C}.\mathrm{L}.).\end{array}$$(42)

Fig. 10. Radius-mass diagram for the 12 EOSs used in this paper to accomplish a good coverage of the range from soft to stiff EOSs, while still being in agreement with the tidal deformability test in the GW170817 LIGO/Virgo event (Abbott et al. 2017, 2018) (black curves). The blue dashed line indicates the mass of PSR J0337+1715, and the red dotted line corresponds to the most massive neutron star used in our combined test: PSR J0348+0432 (Antoniadis et al. 2013). The black curves correspond to the following EOSs (from left to right in their intersection with the red dotted line): WFF1, SLy4, WFF2, AP4, BSk20, ENG, SLy9, AP3, BSk25, BSk21, MPA1, BSk22 (see Lattimer & Prakash 2001 and https://compose.obspm.fr). Our stiffest EOS, BSk22, is also in agreement with the (more model dependent) constraint of Mmax ≲ 2.3 M⊙ by Rezzolla et al. (2018), Shibata et al. (2019). We have included EOS H4 (grey curve), which is disfavoured by GW170817, and therefore has not been used in Figs. 11 and 12. All these EOSs also agree with the latest constraints from NICER (Miller et al. 2019).

Consequently, using Eq. (41), while keeping in mind that according to Eq. (38) ζ ≥ 0, one finds

$$\begin{array}{c}\hfill {\omega }_{\mathrm{BD}}>140000(95\%\mathrm{C}.\mathrm{L}.).\end{array}$$(43)

This limit is more than a factor of three larger, that is, more constraining, than the Cassini limit (Will 2018a). When using EOS H4, which is already disfavoured by the GW170817 LIGO/Virgo event, we find ω_BD >  130 000, which is only marginally weaker than the above limit. Just to illustrate the EOS dependence of the limit on JFBD gravity, for the soft EOS WFF1 (outer/left in Fig. 10), the lower limit for ω_BD increases to 180 000.

While a stiffer EOS gives a more conservative limit for JFBD gravity, such a general statement is no longer true for the whole α₀–β₀ parameter space of T₁(α₀, β₀). In particular for certain negative values of β₀, a softer EOS can be more conservative. For low and medium mass neutron stars, like PSR J0337+1715, the β₀ range where that is the case is rather small (see Fig. 11). For high mass neutron stars the situation is quite different. A soft EOS that has a maximum mass close to the mass of the neutron star leads to considerably weaker limits for all β₀ ≲ −2 (Shibata et al. 2014; Shao et al. 2017). This is of particular importance for constraints from pulsars like PSR J0348+0432 (Antoniadis et al. 2013) (see Fig. 10). Hence for our combined constraints on the parameter space of T₁(α₀, β₀) theories we used a set of EOSs that provide a good coverage of the range from soft to stiff. Furthermore, if for a given point (α₀, β₀), which corresponds to a specific gravity theory, there is a single EOS from our set with which all pulsar constraints are fulfilled then this point in the theory space is not excluded. For our joint analysis we have used the UFF results from this paper in combination with the dipolar radiation tests of PSRs J1012+5307 (Desvignes et al. 2016; Antoniadis et al. 2016), J1141−6545 (Bhat et al. 2008), J1738+0333 (Freire et al. 2012b; Zhu et al. 2019), J1909−3744 (Desvignes et al. 2016; Arzoumanian et al. 2018), and J2222−0137 (Cognard et al. 2017). Our results are shown in Fig. 12.

Fig. 11. PSR J0337+1715 constraints on the α0–β0 space of scalar-tensor theories (Damour & Esposito-Farèse 1993, 1996) from Eq. (17): the area under the curve is still allowed by experiments. Two different neutron-star equations of state are used: a soft one, WFF1 (red curve), and a stiff one, BSk22 (blue curve). The two solid lines use the SEP constraint of this paper. The grey curves show the 2σ-limits from Solar-system experiment: Cassini (solid), LLR (dashed), MESSENGER (dotted). JFBD gravity corresponds to β0 = 0 (thin vertical line).

Fig. 12. Combined EOS-agnostic pulsar constraints on the α0–β0 space of scalar-tensor theories (Damour & Esposito-Farèse 1993, 1996) from Eq. (17): the area under the curve is still allowed by experiments. The blue curve is the result of a combination of dipolar radiation tests from six pulsar-WD binary systems (see text for details). The red curve indicates the improvement when the constraints from this paper are added to the dipolar radiation tests. The grey curves show the 2σ-limits from Solar-system experiment: Cassini (solid), LLR (dashed), MESSENGER (dotted). JFBD gravity corresponds to β0 = 0 (thin vertical line).

6. Conclusions

We described in this paper a test of the universality of free fall (UFF) with the pulsar in a triple star system, PSR J0337+1715. The result we obtain for the UFF violation parameter for the MSP is Δ = ( + 0.5 ± 1.8) × 10⁻⁶ (95% C.L.), which can be stated as a limit, |Δ| < 2.05 × 10⁻⁶ (also 95% C.L.). This represents 30% improvement over the previous test using the same pulsar (Archibald et al. 2018). Interestingly, although we obtain a similar value for Δ, the nature of our limit is different: the uncertainty reported in this work is statistical, while the result of Archibald et al. (2018) is largely made of a systematic uncertainty which we did not find necessary in the present analysis. This particular difference makes our two limits difficult to compare in absence of a physically motivated model for the systematic bias, but should also provide an independent verification of the solidity of the result.

Furthermore, in a generic approach we also provide limits for three post-Newtonian strong-field parameters of the three-body interaction, and discuss in detail the relevance of these limits. In view of other binary pulsar limits, it seems that these limits might be of interest only in very specific situations.

As for Archibald et al. (2018), our results are fully consistent with the predictions of GR. This limit strongly constrains SEP violation and any alternative theories of gravity that predict a violation of the universality of free fall for self-gravitating masses (GWEP), particularly for neutron stars with masses similar to that of PSR J0337+1715. In this paper we explicitly calculate these constraints for a wide class of gravity theories, and as a part of this derive EOS-independent constraints on the parameter space of quadratic mono-scalar-tensor gravity. Specifically, for the coupling parameter of Jordan-Fierz-Brans-Dicke (JFBD) gravity we find ω_BD >  140 000, which is the so far the tightest limit for this scalar-tensor theory. We also present new constraints for Damour–Esposito-Farèse (DEF) gravity, a quadratic extension of JFBD gravity. We combine our limit with limits from binary pulsar experiments while accounting for uncertainties in our knowledge of the equation of state (EOS) of neutron-star matter.

In what remains of the paper, we make a detailed comparison of this experiment with the best previous constraints on GWEP/SEP violation. In Sect. 6.1 we make a more detailed comparison with the experiment by Archibald et al. (2018). In Sect. 6.2, we compare our experiment to radiative experiments from binary pulsars, which have also produced strong and complementary constraints on GWEP violation via their strong constraints on the emission of dipolar gravitational waves. Furthermore, we compare the present limit with potential future limits on dipolar radiation from binary neutron-star and neutron star-black hole mergers.

In all of these experiments, no GWEP violation can be detected; gravity behaves, to within observable precision, as described by GR, which is conjectured to be the only viable theory which fully embodies the SEP.

6.1. Comparison with previous work on J0337+1715

This work distinguishes itself from the Archibald et al. (2018) on the following points:

1. Independent data set;

2. Independent timing model including additional effects;

3. Statistics-limited versus systematics-limited accuracy;

4. Tension in the mass measurements and the first measurement of Ω_O;

5. Generic test of those strong-field post-Newtonian parameters, which are a-priori unconstrained, even within a broad class of scalar-tensor theories.

6. EOS-agnostic constraints on DEF gravity, while accounting for the latest observational constraints on the range of EOSs.

Point (1) benefits from the well sampled timing data acquired by the Nançay radio telescope alone (Sect. 2). All the observations used here were conducted within the same frequency range (1.2−1.7 GHz) and with the same environmental conditions since Nançay is a meridian Kraus design telescope. The NUPPI instrumentation is also routinely used for long term high-precision timing providing excellent and stable results. The instrumentation did not change since its installation in 2011 and there is no need for any time jump in the whole dataset since 2011.

Point (2) makes use of NUTIMO (Sect. 3) which has the specificity of allowing for a fully self-consistent treatment of astrometric parameters via the binding with TEMPO2 and the inclusion of Kopeikin and Shklovskii delays in the model. The model in Archibald et al. (2018) did not include these delays and used a local linear approximation for astrometric corrections. The argument in favour of such proxy was that any systematic effect caused by these approximations should not affect the main SEP signature which has a different frequency. However, we observe that the astrometry then found differs significantly from prior knowledge and in particular Gaia observations which led Archibald et al. (2018) to acknowledge that the resulting astrometry should not be used for other applications. In addition, NUTIMO also fits consistently for DM and DM variations and allows to check for local epoch DM variations (DMX parameters in TEMPO2) which revealed no fluctuations over time. Note that Archibald et al. (2018) did fit DM over 1 year time intervals and marginalised over these parameters using the solution of a least-square fit. We also included the aberration delay that neither Ransom et al. (2014) nor Archibald et al. (2018) mention. This delay has a very small amplitude (sum of ∼30 ns sinusoid at the outer period and a ∼0.1 μs sinusoid at the inner period) and therefore can easily be absorbed by other parameter in a fit, but still creates a signal of magnitude larger than the expected SEP sensitivity. In NUTIMO, potential systematic effects that may not be accounted for by the model are absorbed in a re-scaling of the error bars of the times of arrival via the EFAC parameter which ensures a reduced χ² equal to unity. This in turns conservatively increases uncertainties on the posterior parameters.

Point (3) arises from the fact that, in Archibald et al. (2018), most of the total reported uncertainty of Δ (±0.74 × 10⁻⁶ at 68% C.L.) is associated with systematic uncertainties while in this work our uncertainty is mostly statistical. We account for unmodelled systematic effects, mostly a red-noise component, via the EFAC parameter which is responsible for a modest and conservative widening of ∼8% of the uncertainties. The statistical uncertainty of Archibald et al. (2018) is estimated using MCMC sampling similarly as we do in this work and results in Δ = (− 1.1 ± 0.2) × 10⁻⁶ (68% C.L.). Taken alone, this would signify a 5-sigma SEP violation. However, the authors argue that most of the uncertainty comes from unaccounted systematic effects which could generate a signal at the signature frequency of an SEP violation. In other words, it is claimed that the accuracy is systematics-limited while in this work we are statistics-limited. The physical mechanism of the systematics being unknown, Archibald et al. (2018) propose to model systematics using an empirical stochastic model where the extraneous signal is a weighted sum of sine and cosine functions at frequencies kf_i + lf_o (k, l being small integers) whose weights are drawn from a single Gaussian distribution for each particular realisation. In order to sample the distribution of Δ caused by different realisations of the (stochastic) systematics, Archibald et al. (2018) bootstrapped many sets of synthetic data from the model, re-fitted the orbital model, and thus obtained a value of Δ for each synthetic dataset. This estimate heavily depends on the modelling choices for which no physical justification is currently available. It also seems unlikely that systematics should occur at frequencies kf_i + lf_o if the physical mechanism is unrelated to orbital motion (as an SEP violation would be), unless the signal at these frequencies is the tail (in Fourier space) of a systematic signal which peaks at a different frequency, but should then be seen in a periodogram such as Fig. 4. In this respect, Fig. 4 does not suggest that we should consider a systematics-limited regime in the present work. Although some systematics are present as red noise (see Sect. 4.2), there is no sign of an additional signal around the signature frequency.

Point (4) is about comparing the parameter set of Table 1 with the results of Archibald et al. (2018). A direct comparison is not straightforward because of (i) slightly different definition of the parameters (see Sect. 3.1) and (ii) the fact that most of the parameters are not constants of motion but are defined either at the reference time T_ref or T_pos. To minimise the span of numerical computations we do not use the same reference time as Archibald et al. (2018). However, one can compare masses which are constants of motion (and, in the same way, Δ). The values reported in Archibald et al. (2018), m_p = 1.4359(3) M_⊙, m_i = 0.19730(4) M_⊙, m_o = 0.40962(9) M_⊙¹¹, whose statistical 68% confidence intervals are about 5 times better than ours – similar to the ratio between the uncertainties on Δ – are in tension with the values we report in Table 1, with Δm_p ≃ 2.9σ, Δm_i ≃ 2.6σ, Δm_o ≃ 2.4σ where σ is half of the 68% confidence interval reported in this paper, these differences are much more significant than for Δ. Due to the very large correlation between Δ and the orbital parameters on which depend the masses (period and semi-major axis), the systematic uncertainty estimated in Archibald et al. (2018) for Δ should be similar for the masses (but not reported) and partly release the tension.

In addition, we report the first measurement of the outer longitude of ascending node, Ω_O, which was deemed unconstrained in Archibald et al. (2018) although the dispersion of the fit residuals in Archibald et al. (2018) is smaller than ours. We speculate that the absence of Kopeikin delay in their analysis prevented that measurement. However, this parameter is uncorrelated with Δ and should therefore not affect the SEP test, but its absence should bias astrometric parameters.

Point (5) is related to our implementation of the first post-Newtonian equations of motion, derived from the modified Einstein-Infeld-Hoffmann Lagrangian for strongly self-gravitating masses (Appendix A). In this we use two different approaches, a (mostly) generic one where three of the 12 1PN strong-field parameters are unconstrained by Solar System experiments, and a second approach where, under additional assumptions, all the strong-field 1PN parameters are tightly constrained by adopting binary pulsar constraints for the neutron star sensitivity and its derivative. A detailed motivation for the two different approaches is given in Sect. 5. In the first approach we find generic limits for the remaining three 1PN strong-field parameters, which however are not very tight, and therefore generally not of particular interest. Archibald et al. (2018) do not provide an equally generic analysis as done in our first approach.

Point (6) refers to the constraints of quadratic mono-scalar-tensor gravity, where Archibald et al. (2018) have used a single (outdated) EOS. In our combined tests we have fully accounted for our imperfect knowledge of the EOS of neutron-star matter, and used a set of modern EOSs that covers the range from soft to stiff EOSs. A reason for that is the fact that the most conservative pulsar limits do not always come from the stiffest EOS. Our set of EOSs is in agreement with the latest constraints from LIGO/Virgo and NICER, and can account for the largest neutron star masses measured to date (Antoniadis et al. 2013; Cromartie et al. 2019). We would like to point out, that the recent limits of Capano et al. (2020) exclude some of the stiffer EOSs used in our analysis, which consequently leads to even more stringent constraints than the ones shown in Fig. 12, in particular for β₀ ≳ −1. In view of this, our limits can be considered as conservative.

6.2. Comparison with radiative tests

In GR, the lowest source multipole moment that generates gravitational waves is the quadrupole moment (Thorne 1980). In alternatives to GR, however, one finds lower multipoles, where for the dynamics of a binary system, the dipole moment is the most important one. The occurrence of these lower multipoles is closely related to a violation of the SEP (see e.g. Will 2018a, for a discussion). In scalar-tensor theories, for instance, an asymmetry in the sensitivities s_a in a binary system gives rise to scalar dipolar radiation (see Eq. (35)). While a difference in sensitivity is also the reason for a violation of GWEP, where masses with different compactness are falling differently in an external gravitational field (see Eq. (21)). In a sense, the UFF experiment with PSR J0337+1715 and constraints on dipolar radiation damping with binary pulsars are exploring two different sides of the same coin. The limits in Fig. 12 show that currently for a large part of the parameter space, the test with PSR J0337+1715 is more constraining than dipolar radiation tests from binary pulsars. For sufficiently, negative β₀, however, gravitational wave tests with binary pulsars become more constraining, in particular for small α₀. We have a more detailed discussion on this further below.

Gravitational wave observation of a double neutron-star merger can also be used to constrain the emission of dipolar gravitational waves, as has been done for the first LIGO/Virgo binary neutron-star merger GW170817 (Abbott et al. 2019). Limits on scalar-tensor theories as discussed here, from LIGO/Virgo observations, however, are not expected to be competitive with Solar System and pulsar experiments for most of the parameter space (see Shao et al. 2017). Future ground based gravitational wave detectors have the potential to improve on limits presented here, in particular in a β₀ range which is difficult to constrain with pulsar experiments (see Fig. 13). Future gravitational wave observations of mixed (black hole + neutron star) mergers, in particular the combination of multiple events or the combination of ground and space based gravitational-wave observatories promise significant improvements (see e.g. Carson et al. 2020).

Fig. 13. Comparison of the PSR J0337+1715 constraints of this paper with expected constraints from future gravitational wave observatories of a single binary neutron-star merger: Cosmic Explorer (CE; blue dotted curve) and Einstein Telescope (ET; blue dashed curve). CE and ET curves are taken from Fig. 9 in Shao et al. (2017). Like for the CE and ET curves, the PSR J0337+1715 curve is also based on EOS AP4. Solar System constraints (grey) are as in Fig. 11. JFBD gravity corresponds to β0 = 0 (thin vertical line).

As a final comment, there is an important difference between the UFF test conducted with PSR J0337+1715 and dipolar radiation tests. As discussed in the previous subsection, in the regime of spontaneous scalarisation, the neutron star charge can become (almost) independent of the parameter ζ in the sense that $\zeta s_{\text{p}}^2~1$ remains practically fixed for ζ → 0. In such a situation, the effective gravitational constant in the interaction between a neutron star and a white dwarf becomes indistinguishable from G_N and the test with J0337+1715 becomes practically insensitive to such deviations from GR. Dipolar radiation test with pulsar-white dwarf systems, in contrast, are extremely constraining with respect to such scalarisation phenomena, as can be seen from Eq. (35). More generally, in situations where only the strong field of a neutron star can source additional (long-range) gravitational fields that lead to deviations from GR, the UFF test with PSR J0337+1715 cannot place any constraints, in contrast to radiative tests. Hence, both types of tests are complementary and valuable. Binary pulsar tests have already tightly constrained the occurrence of spontaneous scalarisation in neutron stars. However, depending on the EOS and mass of the neutron star, spontaneous scalarisation is not yet fully ruled out by such experiments (Shao et al. 2017).

Supplementary Figures

A high-resolution version of Fig. 2 Access here

Additional figure: chain plot Access here

Acknowledgments

G. Voisin acknowledges support of the European Research Council, under the European Unions Horizon 2020 research and innovation programme (grant agreement No.715051; Spiders). GV would like to thank F. Mottez and R. P. Breton for their support and helpful discussions during this project. GV also thanks A. Archibald for valuable discussions during the very early stages of this work. 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. We acknowledge financial support from the Action Fédératrice PhyFOG funded by Paris Observatory and from the ‘Programme National Gravitation, Références, Astronomie, Métrologie (PNGRAM)’ funded by CNRS/INSU and CNES, France. GD, MK and NW gratefully acknowledge support from European Research Council (ERC) Synergy Grant ‘BlackHoleCam’ Grant Agreement Number 610058. This work made use of the Scipy libraries (www.scipy.org). The authors would like to thank Lijing Shao for his valuable comments that helped to improve the manuscript.

References

Appendix A: Strong-field equations of motion

In our timing model, the motion of the three bodies follows the equations of motion derived from the post-Galiean-invariant N-body Lagrangian of the modified Einstein-Infeld-Hoffmann (mEIH) formalism (Will 1993; Damour & Taylor 1992). The mEIH equations of motion describe the first post-Newtonian dynamics of a N-body system which also contains strongly self-gravitating masses, under the assumption that the gravitational interaction is Poincaré invariant. Furthermore, it assumes that there are no ‘asymmetric’ terms in the Lagrangian, which are anyway absent in many well motivated theories of gravity (see the discussion in Nordtvedt 1985; Damour & Taylor 1992). The mEIH formalsim is a generalisation of the parametrised post-Newtonian (PPN) equations of motion for fully conservative theories with ξ = 0, in order to include effects related to the strong internal fields of strongly self-gravitating objects, like neutron stars. The mEIH Lagrangian can be written as (cf. Damour & Esposito-Farèse 1992)

$$\begin{array}{cc}\hfill \mathcal{L}& ={\displaystyle \sum _{a=1}^n(-m_a{c}^2+m_a\frac{v_a^2}{2}+m_a\frac{v_a^4}{8c^2})}\hfill \\ \hfill & +\frac{1}{2}{\displaystyle \sum _{a=1}^n\sum _{b\ne a}^n\{\frac{G_{\mathit{ab}}{m}_a{m}_b}{r_{\mathit{ab}}}[1-\frac{({\mathit{v}}_a·{\mathit{n}}_{\mathit{ab}})({\mathit{v}}_b·{\mathit{n}}_{\mathit{ab}})}{2c^2}}\hfill \\ \hfill & -\frac{7}{2}\frac{{\mathit{v}}_a·{\mathit{v}}_b}{c^2}+\frac{3}{2}(\frac{v_a^2}{c^2}+\frac{v_b^2}{c^2})+{\overline{\gamma }}_{\mathit{ab}}\frac{{({\mathit{v}}_a-{\mathit{v}}_b)}^2}{c^2}]\hfill \\ \hfill & -{\displaystyle \sum _{c\ne a}^n\frac{G_{\mathit{ab}}{G}_{\mathit{ac}}{m}_a{m}_b{m}_c}{c^2{r}_{\mathit{ab}}{r}_{\mathit{ac}}}(1+2{\overline{\beta }}_{\mathit{bc}}^a)}\},\hfill \end{array}$$(A.1)

where m_a are the inertial masses with coordinate positions x_a and coordinate velocities v_a, r_ab = ∥x_a − x_b∥, v_a = ∥v_a∥, and n_ab = (x_b − x_a)/r_ab. The quantities ${\overline{\gamma }}_{\mathit{ab}}={\gamma }_{\mathit{ab}}-1$, ${\overline{\beta }}_{\mathit{bc}}^a={\beta }_{\mathit{bc}}^a-1$ and G_ab are the effective strong-field interaction constants. The unbarred quantities are the strong-field generalisation of the PPN parameters γ_PPN and β_PPN (Eddington parameters). The strong-field parameters satisfy the symmetries G_ab = G_ba (a ≠ b), ${\overline{\gamma }}_{\mathit{ab}}={\overline{\gamma }}_{\mathit{ba}}$ (a ≠ b), and ${\overline{\beta }}_{\mathit{bc}}^a={\overline{\beta }}_{\mathit{cb}}^a$ (a ≠ b, a ≠ c). The body-dependent effective strong-field interaction constants depend on the details of the underlying gravity theory as well as the structure of the individual bodies. Hence, in the most general case of a three-body system one has three different effective gravitational constants G_ab, three different ${\overline{\gamma }}_{\mathit{ab}}$, and nine different ${\overline{\beta }}_{\mathit{cb}}^a$. In GR, due to the fulfilment of SEP and the corresponding effacement of the internal structure (see e.g. Damour 1987), one has G_ab = G_N and ${\overline{\gamma }}_{\mathit{ab}}={\overline{\beta }}_{\mathit{bc}}^a=0$.

One can then use the Euler-Lagrange equations (Will 1993) to derive the equations of motion for each body:

$$\begin{array}{cc}\hfill {\ddot{\mathit{x}}}_a=& {\displaystyle \sum _{b\ne a}\frac{G_{\mathit{ab}}{m}_b}{r_{\mathit{ab}}^2}{\mathit{n}}_{\mathit{ab}}[1-\frac{1}{c^2}(4{\mathit{v}}_a·{\mathit{v}}_b-{\mathit{v}}_a^2-2{\mathit{v}}_b^2}\hfill \\ \hfill & +\frac{3}{2}{({\mathit{v}}_b·{\mathit{n}}_{\mathit{ab}})}^2-{\overline{\gamma }}_{\mathit{ab}}{({\mathit{v}}_a-{\mathit{v}}_b)}^2)]\hfill \\ \hfill & +{\displaystyle \sum _{b\ne a}\frac{G_{\mathit{ab}}{m}_b}{r_{\mathit{ab}}^2{c}^2}({\mathit{v}}_b-{\mathit{v}}_a)[{\mathit{n}}_{\mathit{ab}}·(4{\mathit{v}}_a-3{\mathit{v}}_b-2{\overline{\gamma }}_{\mathit{ab}}({\mathit{v}}_b-{\mathit{v}}_a))]}\hfill \\ \hfill & +{\displaystyle \sum _{b\ne a}\sum _{c\ne b}\frac{G_{\mathit{ab}}{G}_{\mathit{bc}}{m}_b{m}_c}{r_{\mathit{ab}}{r}_{\mathit{bc}}{c}^2}[\frac{1}{r_{\mathit{bc}}}(\frac{1}{2}({\mathit{n}}_{\mathit{ab}}·{\mathit{n}}_{\mathit{bc}}){\mathit{n}}_{\mathit{ab}}+\frac{7}{2}{\mathit{n}}_{\mathit{bc}})}\hfill \\ \hfill & -\frac{{\mathit{n}}_{\mathit{ab}}}{r_{\mathit{ab}}}+2{\overline{\gamma }}_{\mathit{ab}}\frac{{\mathit{n}}_{\mathit{bc}}}{r_{\mathit{bc}}}-2{\overline{\beta }}_{\mathit{ca}}^b\frac{{\mathit{n}}_{\mathit{ab}}}{r_{\mathit{ab}}}]\hfill \\ \hfill & -{\displaystyle \sum _{b\ne a}\sum _{c\ne a}\frac{G_{\mathit{ab}}{G}_{\mathit{ac}}{m}_b{m}_c}{r_{\mathit{ab}}^2{r}_{\mathit{ac}}{c}^2}{\mathit{n}}_{\mathit{ab}}[4+2{\overline{\gamma }}_{\mathit{ac}}+2{\overline{\beta }}_{\mathit{bc}}^a].}\hfill \end{array}$$(A.2)

In the weak-field limit one can check that this equation does give the PPN equation of motion (e.g. Soffel 1989; Will 1993).

Conserved quantities are key elements to check the numerical implementation and accuracy of the equations of motion. We have used the Hamiltonian (conservation of energy), and the momentum and position of the centre of mass of the system. The last two are also necessary to derive the initial conditions of the system.

The Hamiltonian corresponding to Eq. (A.1) is derived using the Legendre transform $\mathcal{H}=\sum _a{\mathit{v}}_a·\frac{\partial \mathcal{L}}{\partial {\mathit{v}}_a}-\mathcal{L}$,

$$\begin{array}{cc}\hfill \mathcal{H}=& {\displaystyle \sum _a\{m_a{c}^2+\frac{1}{2}{m}_a{v}_a^2+\frac{3}{8}{m}_a\frac{v_a^4}{c^2}}\hfill \\ \hfill & +\frac{1}{2}{\displaystyle \sum _{b\ne a}[-\frac{G_{\mathit{ab}}{m}_a{m}_b}{r_{\mathit{ab}}}(1+\frac{7}{2}\frac{{\mathit{v}}_a·{\mathit{v}}_b}{c^2}}\hfill \\ \hfill & +\frac{1}{2c^2}({\mathit{v}}_a·{\mathit{n}}_{\mathit{ab}})({\mathit{v}}_b·{\mathit{n}}_{\mathit{ab}})-3\frac{v_a^2}{c^2}+{\overline{\gamma }}_{\mathit{ab}}\frac{{({\mathit{v}}_a-{\mathit{v}}_b)}^2}{c^2})\hfill \\ \hfill & +\frac{1}{2}{\displaystyle \sum _{c\ne a}(1+2{\overline{\beta }}_{\mathit{bc}}^a)\frac{G_{\mathit{ab}}{m}_a{m}_b}{r_{\mathit{ab}}}\frac{G_{\mathit{ac}}{m}_c}{r_{\mathit{ac}}{c}^2}}]\}·\hfill \end{array}$$(A.3)

The momentum of the centre of mass is given by the same expression as in GR only with the replacement G → G_ab, $\mathit{P}=\sum _a\frac{\partial \mathcal{L}}{\partial {\mathit{v}}_a}$ and

$$\begin{array}{cc}\hfill \mathit{P}=& {\displaystyle \sum _a[m_a{\mathit{v}}_a(1+\frac{v_a^2}{2c^2}-\frac{1}{2}\sum _{b\ne a}\frac{G_{\mathit{ab}}{m}_b}{c^2{r}_{\mathit{ab}}})}\hfill \\ \hfill & -\frac{1}{2}{\displaystyle \sum _{b\ne a}\frac{G_{\mathit{ab}}{m}_b}{c^2{r}_{\mathit{ab}}}({\mathit{v}}_a·{\mathit{n}}_{\mathit{ab}}){\mathit{n}}_{\mathit{ab}}}].\hfill \end{array}$$(A.4)

The centre-of-mass position X satisfies $(\mathcal{H}/c^2)\frac{\mathrm{d}X}{\mathrm{d}t}=\mathit{P}$ (see e.g. Will 2014b),

$$\begin{array}{c}\hfill \frac{\mathcal{H}}{c^2}\mathit{X}={\displaystyle \sum _a{m}_a{\mathit{x}}_a(1+\frac{v_a^2}{2c^2}-\frac{1}{2}\sum _{b\ne a}\frac{G_{\mathit{ab}}{m}_b}{c^2{r}_{\mathit{ab}}})+\mathcal{O}(c^{-4}).}\end{array}$$(A.5)

All Tables

All Figures

	Fig. 1. Top: template pulse profile of PSR J0337+1715 used for timing. 450 h of observations conducted with the Nançay Radio Telescope were integrated over the frequency range 1230–1742 MHz. Bottom: pulse profile obtained on October 4, 2014. When profiles are ToA uncertainty ranked, this one is at the 10th percentile of the lowest uncertainties, typical of a good observation.
In the text

	Fig. 2. Corner plot of the correlations between the fitted parameters of Table 1 and their marginalised distribution (diagonal histograms) sampled by our MCMC. A high-resolution version is available as online supplementary material.
In the text

Fig. 3. Sketch of the orbits of PSR J0337+1715 (not to scale). Note that the orbits are nearly circular and that the ellipsoidal shape of the orbits on the left-hand sketch arises purely from projection. The osculating orbits of the system can be parametrised by the Keplerian orbital elements of the pulsar within the inner binary and of the inner binary within the outer binary (see text). In particular, ap is the semi-major axis of the pulsar orbit within the inner binary, ab is the semi-major axis of the inner binary within the outer binary and eO its eccentricity. Are also shown the longitude of periastron the outer binary as well as the inner and outer inclinations with respect to the plane of the sky, iI and iO respectively. Note that for simplicity the sketch neglects the difference between the directions of the inner and outer lines of ascending nodes, δΩ, as it is in practice very small.

In the text

Fig. 4. Lomb-Scargle periodogram of the timing residuals for the solution of Table 1. The periodogram is sampled at fm/5 where fm ≃ 2268 days−1 is the inverse of the full time span. From left to right vertical lines show the frequencies of the red-noise component (fRN ≃ 1650 day−1, black), the Earth orbital period (fE, green), the outer-orbit orbital period (fO, orange), the inner-binary orbital period (fI, red) and the SEP violation signature (2fI − fO, purple).

In the text

	Fig. 5. Relative error in energy conservation during the numerical integration of the equations of motion over the 2268 day span of our data. The envelope of the signal oscillates at the outer binary frequency while a zoom would show a fast oscillation at the inner binary orbital frequency.
In the text

Fig. 6. Three components of the centre-of-mass position variation during the integration of the equations of motion over the 2268 day span of our data. A quadratic component has been fitted out (see text and legend) which leaves only the oscillatory components. The main visible oscillation is at the frequency of the outer binary while a zoom would show a fast oscillation at the inner binary orbital frequency.

In the text

	Fig. 7. Left-hand side: marginalised posterior probability distribution of the SEP violation parameter Δ sampled by MCMC and normal law with the same mean and standard deviation. The upper axis gives the mean value and the boundaries of the 95% confidence region. Right-hand side: distribution of distance to GR derived from the left-hand-side distribution.
In the text

	Fig. 8. Residuals of the times of arrival using the mean parameters reported in Table 1.
In the text

Fig. 9. Residuals of the times of arrival corresponding to the parameters of Table 1 plotted versus the inner-binary orbital phase (top), the outer-binary orbital phase (middle) and the Earth orbital phase (bottom). The red vertical bars of the first two plots show the phase where the pulsar (resp. the inner binary) is closer to the Solar System. The red vertical bar on the bottom plot, shows the Earth orbital phase where the line of sight to the pulsar passes closest to the Sun. We have removed every ToA when the pulsar is less than 5 deg from the Sun, which explains the gap around that particular phase.

In the text

Fig. 10. Radius-mass diagram for the 12 EOSs used in this paper to accomplish a good coverage of the range from soft to stiff EOSs, while still being in agreement with the tidal deformability test in the GW170817 LIGO/Virgo event (Abbott et al. 2017, 2018) (black curves). The blue dashed line indicates the mass of PSR J0337+1715, and the red dotted line corresponds to the most massive neutron star used in our combined test: PSR J0348+0432 (Antoniadis et al. 2013). The black curves correspond to the following EOSs (from left to right in their intersection with the red dotted line): WFF1, SLy4, WFF2, AP4, BSk20, ENG, SLy9, AP3, BSk25, BSk21, MPA1, BSk22 (see Lattimer & Prakash 2001 and https://compose.obspm.fr). Our stiffest EOS, BSk22, is also in agreement with the (more model dependent) constraint of Mmax ≲ 2.3 M⊙ by Rezzolla et al. (2018), Shibata et al. (2019). We have included EOS H4 (grey curve), which is disfavoured by GW170817, and therefore has not been used in Figs. 11 and 12. All these EOSs also agree with the latest constraints from NICER (Miller et al. 2019).

In the text

Fig. 11. PSR J0337+1715 constraints on the α0–β0 space of scalar-tensor theories (Damour & Esposito-Farèse 1993, 1996) from Eq. (17): the area under the curve is still allowed by experiments. Two different neutron-star equations of state are used: a soft one, WFF1 (red curve), and a stiff one, BSk22 (blue curve). The two solid lines use the SEP constraint of this paper. The grey curves show the 2σ-limits from Solar-system experiment: Cassini (solid), LLR (dashed), MESSENGER (dotted). JFBD gravity corresponds to β0 = 0 (thin vertical line).

In the text

Fig. 12. Combined EOS-agnostic pulsar constraints on the α0–β0 space of scalar-tensor theories (Damour & Esposito-Farèse 1993, 1996) from Eq. (17): the area under the curve is still allowed by experiments. The blue curve is the result of a combination of dipolar radiation tests from six pulsar-WD binary systems (see text for details). The red curve indicates the improvement when the constraints from this paper are added to the dipolar radiation tests. The grey curves show the 2σ-limits from Solar-system experiment: Cassini (solid), LLR (dashed), MESSENGER (dotted). JFBD gravity corresponds to β0 = 0 (thin vertical line).

In the text

Fig. 13. Comparison of the PSR J0337+1715 constraints of this paper with expected constraints from future gravitational wave observatories of a single binary neutron-star merger: Cosmic Explorer (CE; blue dotted curve) and Einstein Telescope (ET; blue dashed curve). CE and ET curves are taken from Fig. 9 in Shao et al. (2017). Like for the CE and ET curves, the PSR J0337+1715 curve is also based on EOS AP4. Solar System constraints (grey) are as in Fig. 11. JFBD gravity corresponds to β0 = 0 (thin vertical line).

In the text