Guillaume Voisin
email: guillaume.voisin@manchester.ac.uk

Ismaël Cognard

Paulo C C Freire

Norbert Wex

Lucas Guillemot

Grégory Desvignes

Michael S Kramer

Gilles Theureau

An improved test of the strong equivalence principle with the pulsar in a triple star system

Keywords: Gravitation -(Stars, ) pulsars, individual PSR J0337+1715 -Stars, neutron -Radio continuum, stars

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Newton | Philosophiae Naturalis Principia Mathematica[END_REF]. 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 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 2018b). 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 2018b).

This generalisation of the relativity principle to reference frames in free fall guided Einstein towards general relativity (GR, [START_REF] Einstein | Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften[END_REF]. 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 selfgravitational 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 2018b) 1 .

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;Shao et al. 2013;Shao & Wex 2013). In what follows, we focus on the gravitational version of the WEP (GWEP, Will 2018b), 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 ẍa = -G ab m b r ab r ab 3 + O(c -2 ) ,

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:

G ab = m P m a m A m b G N , (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 semiconservative metric theories of gravity that have a conservation of momentum one has only a single gravitational mass m G ≡ m P = m A (Will 2018b). For the remainder of the paper we assume that momentum is conserved in the gravitational interaction, and therefore G ab = G ba . 2 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

∆ E -∆ M m G m E - m G m M = (-3.0 ± 5.0) × 10 -14
(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 selfgravitating, however, this is especially true for the two 'proof masses', the Earth and the Moon: for the Earth ε grav,E ≡ E grav,E /m E c 2 = -4.6 × 10 -10 (here E grav,E is the Newtonian gravitational binding energy of the Earth), for the Moon ε grav,M = -0.2 × 10 -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 G. Voisin et al.: An improved test of the strong equivalence principle with the pulsar in a triple star system small that it can only have a very small effect on ∆ a . Within the parametrised post-Newtonian (PPN) formalism for metric theories of gravity

∆ a = η ε grav,a , (4)
where η is the so called Nordtvedt parameter, a combination of several PPN parameters (see Will 2018b, for details). The current limit on the Nordtvedt parameter from LLR is (-0.2±1.1)×10 -4 . 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, [START_REF] Genova | [END_REF] have derived η = (-6.6 ± 7.2) × 10 -5 .

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: 3 Stairs et al. (2005) derive |∆| < 5.6 × 10 -3 , and Gonzalez et al. (2011) derive |∆| < 4.6 × 10 -3 (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 [START_REF] Wex | [END_REF].

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 [START_REF] Damour | Physics of Relativistic Objects in Compact Binaries: From Birth to Coalescence[END_REF].

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 -3 (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 strongfield 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 -10 m s -2 in the solar vicinity. In the case of the LLR test, the polarising gravitational field (that of the Sun) is much stronger (6 × 10 -3 m s -2 ), 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 -3 m s -2 ), 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 driftscan 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.4M 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 [START_REF] Ransom | [END_REF]. 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 [START_REF] Tauris | [END_REF].

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 [START_REF] Ransom | [END_REF]. More importantly, the GWEP test was eventually carried out for PSR J0337+1715 by Archibald et al. (2018), yielding |∆| ≤ 2.6 × 10 -6 (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 monoscalar-tensor theories described by Damour & Esposito-Farèse (1992, 1993, henceforth DEF gravity). The constraints on the weak-field coupling parameter for these theories (α 0 ) derived from PSR J0337+1715 significantly improve upon all previous experiments for most of their β 0 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 2018b); 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 Section 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 Section 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 Section 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 Section 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 2018b). We also derive new constraints on a sub-class of those theories, the Damour-Esposito Farèse (DEF) theories (Damour & Esposito-Farèse 1992;Damour & Esposito-Farèse 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 Section 6.

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-meter dish able to conduct ∼1 hour observations on any given source within its declination range each day. The dual linear polarisation signals are sent to the Nançay Ultimate Pulsar Processing Instrument (NUPPI, Desvignes et al. 2011), 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, Du-Plain 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-second sub-integrations. A single standard timing parameter file in tempo format4 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) 5 . A ∼3 Hz pulsed noise diode is fired for ∼10 seconds 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-second 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 minutes, 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 r.m.s. 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).

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 [START_REF] Ransom | [END_REF] 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 ∼ Gm ac 2 ∼ v 2 c 2 5 × 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 [START_REF] Will | Theory and experiment in gravitational physics[END_REF] and Damour & Taylor (1992) which parametrises almost the entire class of 'fully conservative' Lagrangian-based theories of gravity (without preferred location effects6 ) 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 strongfield 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 [START_REF] Will | Theory and experiment in gravitational physics[END_REF] 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 -4 . 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 Section 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: βp , βpp , βp . Since the limits we find for these parameters in Section 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 7 , solves numerically the 3-body equations of motion at 1PN (see appendix A.1 and particularly equation (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 Roemer, Kopeikin and Shklovskii delays (Shklovskii 1970;Kopeikin 1996), and adds an extra second order correction, that is, L 2 /d 2 , for the latter (the only second order correction that may become important with time, see e.g. [START_REF] Voisin | Theses[END_REF]). We note that Kopeikin and Shklovskii delays were not included in previous works (Archibald et al. 2018;[START_REF] Ransom | [END_REF]. 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 2 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 pre-viously 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 Chapter 5 of [START_REF] Voisin | Theses[END_REF].

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 Figure 5). 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.

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 Chapter 5 of [START_REF] Voisin | Theses[END_REF] for details about the delays). This is a common practice in pulsar timing. In addition, we rescaled 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 Eqn. (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 f and f are connected to the intrinsic parameters f and f by

f = f 1 - v 2 + U i + U o 2c 2 -µ 2 ⊥ d∆T 1 - 3 2 ∆T µ d /c , (5)
f = f -f dµ 2 ⊥ (1 -3∆T µ d )/c , (6)
where

v = 2π a b P O + a p P I , (7)
U i = Gm i a p (1 + m p /m i ) , (8)
U o = Gm o a b (1 + (m p + m i )/m o ) , (9)
∆T = T ref -T pos , (10)
where the symbols correspond to those defined in Table 2. The use of the two re-scaled parameters above instead of the intrinsic Fig. 2. 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, a p is the semi-major axis of the pulsar orbit within the inner binary, a b is the semi-major axis of the inner binary within the outer binary and e O 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, i I and i O 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.

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 Figure 2). 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.

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 Ii 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 G. Voisin et al.: An improved test of the strong equivalence principle with the pulsar in a triple star system 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 8 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 = (H /c 2 )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 ).

v i = v i/b +v b .

Model accuracy

The main signature of a SEP violation in our pulsar-timing experiment would be a residual signal at the frequency 2 f if o (Archibald et al. 2018), where f i,o are the inner and outer orbital frequency respectively (see also Figure 9). 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 if 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 2 f if o as originally pointed out by Archibald et al. (2018). Using the Newtonian-order linear theory of Nordtvedt (1968), one can show that |∆| ∼ 10 -6 translates into a signal amplitude of ∼ 0.1 s via the Roemer 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 Figure 5). 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).

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 [START_REF] Stoer | Introduction to numerical analysis[END_REF] implemented in the Odeint module [START_REF] Ahnert | NUMERICAL ANALY-SIS AND APPLIED MATHEMATICS ICNAAM 2011: International Conference on Numerical Analysis and Applied Mathematics[END_REF]) of the Boost library9 . 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 roundoff errors. In practice we observe round-off errors at the level of 10 -3 ns showing that numerical round-off errors are not an issue given our objective of a 1ns 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 [START_REF] Press | Numerical recipes in FORTRAN[END_REF]) 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 χ 2 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 ∼ 1ns signal, then the χ 2 can be expanded as follow,

χ 2 = N i=1 ∆t 2 i σ 2 i + N i=1 ∆t i δ i σ 2 i + N i=1 δ 2 i σ 2 i , ( 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 χ 2 ∼ N. The second term in equation ( 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δ 2 /σ 2 . Taking δ i ∼ δ = 1 ns we see that the second and third term of equation (11) are respectively 5 × 10 -4 χ 2 and ∼ 2.5 × 10 -7 χ 2 . 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 χ 2 between two increments is smaller than 2.5 × 10 -7 .

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 Section 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 2 /c 2 ∼ 5 × 10 -7 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 accu- rate to first order (1PN) and that systematic 'residuals' of order v 4 /c 4 ∼ 2 × 10 -13 (2PN) are present in the equations themselves. As a result, we see on Figure 3 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 ẍ → ẋ → 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 3 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 Figure 4 would bias the spin frequency by ∼ 10 -14 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 Figure 4. The residuals show that the systematic oscillatory 2PN motion 4.8 mas yr -1 1 mas yr -1 † 1 µ δ -4.4 mas yr -1 0.8 mas yr -1 † 1 V = µ d d 29.7 km s -1 0.9 km s -1 * 2 † : 2× the uncertainty reported in the source. * : 1× the uncertainty reported in the source.

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.

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. 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 1).

Bayesian analysis results

Our

We have also applied a factor of two to the uncertainty on the photometric distance reported in [START_REF] Ransom | [END_REF] 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 [START_REF] Yao | [END_REF], both predict a distance of about 800 pc which is significantly smaller that reported in Tables 12. This indicates that the electron density for Notes. The error bars apply to the digits between parenthesis. Upper-case indices I, O refer to the inner and outer binary respectively while lower-case indices p, i, o refer to the pulsar, inner white dwarf and outer white dwarf respectively. the given Galactic height (z = -690(40) pc) is overestimated. All priors are summarised in Table 1. Let us note that our fit for the radial proper motion µ d is unconstraining as the uncertainties reported in Table 2 match the radial velocity prior of Table 1. 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 affineinvariant 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 [START_REF] Allison | [END_REF]) 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 Figure 5). 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 Section 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), Section 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, βp , βpp , βp , yielding the following 95% C. L. constraints for them:

∆ = 1.1 +4.8 -4.6 × 10 -6 , (13) βp = -2 +13 -14 , (14)
βpp = -0.1 +0.9 -1.0 , (15) βp = +0.8 +7.2 -6.6 .

In regard to binary-pulsar tests (see Section 5), the above results on the three β parameters are unconstraining. We used this prior knowledge to run our main model with βp = βp = βpp = 0 and obtain our primary SEP limit

∆ = (+0.5 ± 1.8) × 10 -6 (95% C.L.), (17)
which translates into |∆| < 2.05 × 10 -6 at 95% C. L. (see Figure 6). The full result of the main model is reported in Table 2. Note that ∼ 8% of the reported uncertainties are due to unaccounted systematics absorbed in the EFAC parameter (see also Section 4.2). The wider uncertainty obtained in the secondary run is due to large correlations with the three additional parameters.

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 Section 3.1) but very large linear correlations, as can be seen in Figure 5, 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 Section 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 χ 2 (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),

rô= m, σ = σ ( ô) σ , ( 18
)
where σ is a standard-deviation estimator, ô is a statistical estimator which here is either the mean m or the standard deviation σ, 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 69408 recorded elements (241 independent sets of 288 walkers). The standard deviation of ô was estimated by applying ô 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 r m and r σ 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.

Analysis of the residuals

We have assessed the robustness of our fit by evaluating the presence of systematic components in the residuals (Figure 7). As it appears from Figure 8, 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.

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, Figure 9. 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(Hobbs et al. , 2012) ) or inaccuracies in the planetary masses used to generate the Solar System ephemeris [START_REF] Champion | [END_REF]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 [START_REF] Champion | [END_REF]. 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 2) in order to accommodate this systematic effect into a reduced χ 2 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 -1 and its harmonics on the periodogram of Figure 9 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 ∼ 1h windows separated by an integer number of sidereal day. Due to the proximity of the inner orbital frequency with the observing frequency of ∼ 1day -1 one might be concerned with a potential bias. However we have checked that the Fourier transform of a comb of 1h 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.

Theory dependent tests

The parametrised post-Newtonian (PPN) formalism (see e.g. Will 2018b), 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 2018a). 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. [START_REF] Damour | Physics of Relativistic Objects in Compact Binaries: From Birth to Coalescence[END_REF]). Scalartensor theories of gravity have turned out to be particularly useful for this purpose. Apart from being well motivated and well studied alternatives to GR [START_REF] Fujii | The Scalar-Tensor Theory of Gravitation[END_REF], 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;Damour & Esposito-Farèse 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 2018b). Quite a number of well known scalartensor theories belong to this class, like Jordan-Fierz-Brans-Dicke (JFBD) gravity [START_REF] Jordan | Schwerkraft und Weltall[END_REF][START_REF] Fierz | [END_REF]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 Section 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. 2.

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

S = 1 16πG 0 √ -g d 4 x φR - ω(φ) φ ∂ µ φ∂ µ φ -U(φ) +S mat ψ A mat , g µν , ( 19
)
where G 0 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

G N = G 0 φ 0 (1 -ζ) , (20)
where φ 0 denotes the cosmological background field and ζ ≡ 1/(2ω(φ 0 ) + 4). S mat is the action of the matter fields ψ A mat , 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 (φ) 1/a 2 b . 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

G ab = G N 1 -2ζ(s a + s b -2s a s b ) , (21)
where the sensitivity

s a ≡ d ln m a (φ) d ln φ φ 0 (22)
accounts for the dependence of each body on a change in the ambient scalar field, while the number of baryons remains fixed (Will 2018b). 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

∆ ab = -2ζ(s a + s b -2s a s b ). ( 23
)
Because of the product s a s 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 2018b). For weakly self-gravitating bodies the sensitivity s a is simply related to the fractional gravitational binding energy ε grav,a via

s a = -(1 + 2λ) ε grav,a + O ε 2 grav,a , (24)
where

λ ≡ φ 0 ω (φ 0 )ζ 2 (1 -ζ) -1 .
The two Eddington parameters of the PPN formalism are given by

γ PPN = 1 -2ζ , β PPN = 1 + ζλ , (25)
and the Nordtvedt parameter of equation ( 4) reads

η = 4β PPN -γ PPN -3 = 2ζ(1 + 2λ) . ( 26
)
For the inner and outer white dwarf we have ε grav,i -1.8 × 10 -5 and ε grav,o -4.8 × 10 -5 , 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 (2) leads to a direct constraint for

ζ s p - ∆ 2 = (-0.2 ± 0.9) × 10 -6 (95% C.L.) , (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.

First post-Newtonian contributions

At the first post-Newtonian (1PN) level (order v 2 /c 2 ), the PPN parameters βPPN ≡ β PPN -1 and γPPN ≡ γ PPN -1 need to be replaced by the body-dependent quantities γab and βa bc (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 [START_REF] Will | Theory and experiment in gravitational physics[END_REF] and Damour & Esposito-Farèse (1992). The latter uses the so-called Einstein-frame representation and gives these terms for multiscalar-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, γio γPPN ,

(28) γpi γpo γpb (s b = 0) ≡ γp , (29) βi oo βo ii βPPN , (30) βi
po βo pi βa pc (s a = s c = 0) ≡ βp , (31) βi
pp βo pp βa pp (s a = 0) ≡ βpp , (32) βp
ii βp io βp oo βp bc (s b = s c = 0) ≡ βp , (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: γab = γba and βa bc = βa cb . 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 βPPN and γPPN from Solar System tests ∼ 10 -5 (Will 2018b), directly put tight constraints on two of the six 1PN parameters. Furthermore, γp only depends on terms proportional to ζ and ζ s p , the first being constrained to ∼ 10 -5 by Cassini (Bertotti et al. 2003) and the latter to ∼ 10 -6 already by the Newtonian-level dynamics of the triple system (cf. equation (27); see also limit ( 13)). Consequently, without loss of generality, γab , βi oo , and βo ii 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 βp , βpp , and β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 (∆, βp , βpp , βp ), which is called secondary model in Section 3 and 4. Our analysis based on this model leads to the generic limits ( 13) and ( 14) -( 16).

In our generic approach, the parameters βp , β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 β parameters have terms, which are proportional to λζ, ζ 2 s 2 p , λζ s p , λζ s 2 p , and ζ 2 s p , where

s a ≡ d 2 ln m a (φ) d(ln φ) 2 φ 0 (34)
(see e.g. Will 2018a, for details). Solar System constraints on β PPN and η put tight constraints (∼ 10 -5 ) on λζ, and ζ 2 s 2 p is constrained to ∼ 10 -12 because of equation ( 27). However, the quantities λ, s p , and s p are unconstrained by above considerations. Consequently, λζ s p , λζ s 2 p , and ζ 2 s p are a-priori unconstrained. Below, in a more theory specific context, we discuss a situation in DEF gravity where ζ s 2 p 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 p can assume very large values. As a consequence, the β 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 βpp / βPPN ∼ ζ -1 while β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 β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 equation ( 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 ζ s 2 p , 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 Ṗdipole

b - G N c 3 16π 2 P b m p m c m p + m c 1 + e 2 /2 (1 -e 2 ) 5/2 ζ s 2 p ( 35)
(see e.g. equation (12.32) in Will 2018b), 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 ).

ζ s 2 p = (-0.5 ± 1.7) × 10 -6 (95% C.L.) . (36)
As a result, βpp can, in general, also be assumed to be small in the PSR J0337+1715 system. Imposing a generic constraint on ζ 2 s p , and therefore on βp , is somewhat less direct. s p enters, for instance, the precession of periastron, ω, which is particularly well tested -in combination with other post-Keplerian parameters -in eccentric shortorbital-period binary pulsars [START_REF] Wex | [END_REF]Will 2018b). However, only the so called Double Pulsar allows for a generic constraint on deviations from GR in ω, which is of the order of 10 -3 (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 ζ 2 s p , and therefore βp can generally be assumed to be small as well. Moreover, given that the Cassini experiment already imposes ζ 2 10 -10 , s p would have to assume quite extreme values to lead to a significant β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 β parameters are tightly constrained by a combination of equation ( 27) with constraints from binary pulsars experiments. The limits in equations ( 14) to ( 16) are therefore generally not of particular interest, at least for constraining the class of scalar-tensor theories considered in this section.

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

In order to explicitly calculate s p and s p 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 1 (α 0 , β 0 ) 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

ω(φ) = 1 2       1 α 2 0 -β 0 ln φ -3       , (37)
where φ 0 ≡ 1, without loss of generality (Will 2018b). Furthermore, one finds

ζ = α 2 0 1 + α 2 0 , λ = β 0 2(1 + α 2 0 ) . ( 38
)
The tight constraints on ζ from the Cassini mission imply that α 2 0 10 -5 . Furthermore, from equation ( 26) one then finds for the Nordtvedt parameter

η 2α 2 0 (1 + β 0 ) . (39)
The sensitivity s a of a weakly self-gravitating body, like the WD companions to J0337+1715, can be calculated according to equation ( 24):

s a -(1 + β 0 ) ε grav,a . (40)
For neutron stars, the absolute value of the sensitivity can become very large if β 0 -4.5, a fact first discovered within T 1 (α 0 , β 0 ) gravity theories by Damour & Esposito-Farèse (1993), and generally referred to as 'spontaneous scalarisation'. As a result, even for arbitrarily small α 0 , the quantity α 0 s a ζ 1/2 s a remains at order unity. 10 A special case of T 1 (α 0 , β 0 ) is JFBD gravity, for which β 0 = 0. In that case λ = 0, and the coupling function is a constant:

ω(φ) = 1 2α 2 0 - 3 2 = 1 2ζ -2 ≡ ω BD > - 3 2 , (41)
which is called the Brans-Dicke parameter. For ω BD → ∞, that is, α 0 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 2 , this value is practically independent of ζ (Shibata et al. 2014;Shao et al. 2017). Hence, equation ( 27) can directly be converted into limits on the coupling parameter:

ζ = (-1 ± 6) × 10 -6 (95% C.L.) .
(42) 10 The effective scalar coupling α A used by Damour & Esposito-Farèse (1993) is linked to the sensitivity as defined here via α A = α 0 (1 -2s A ). Furthermore, α 0 ζ 1/2 for small α 0 .

Consequently, using equation ( 41), while keeping in mind that according to equation (38) ζ ≥ 0, one finds ω BD > 140 000 (95% C.L.) .

This limit is more than a factor of three larger, that is, more constraining, than the Cassini limit (Will 2018b). 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. (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).

While a stiffer EOS gives a more conservative limit for JFBD gravity, such a general statement is no longer true for the whole α 0 -β 0 parameter space of T 1 (α 0 , β 0 ). In particular for certain negative values of β 0 , a softer EOS can be more conservative. For low and medium mass neutron stars, like PSR J0337+1715, the β 0 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 β 0 -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 1 (α 0 , β 0 ) 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 (α 0 , β 0 ), 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 scalartensor theories (Damour & Esposito-Farèse 1993;Damour & Esposito-Farèse 1996) from equation ( 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).

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 -6 (95% C.L.), which can be stated as a limit, |∆| < 2.05 × 10 -6 (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 This 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 Fig. 12. Combined EOS-agnostic pulsar constraints on the α 0 -β 0 space of scalar-tensor theories (Damour & Esposito-Farèse 1993;Damour & Esposito-Farèse 1996) from equation ( 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).

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 neutronstar 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 Section 6.1 we make a more detailed comparison with the experiment by Archibald et al. (2018). In Section 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 starblack 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.

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 (Section 2). All the observations used here were conducted within the same frequency range (1.2-1.7GHz) 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 (Section 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 [START_REF] Ransom | [END_REF] 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 χ 2 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 -6 at 68% CL) 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 -6 (68% CL). 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 k f i + l f 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 k f i + l f 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 Figure 9. In this respect, Figure 9 does not suggest that we should consider a systematicslimited regime in the present work. Although some systematics are present as red noise (see Section 4.2), there is no sign of an additional signal around the signature frequency.

Point 4) is about comparing the parameter set of Table 2 with the results of Archibald et al. (2018). A direct comparison is not straightforward because of i) slightly different definition of the parameters (see section 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)M11 , 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 2, 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 selfgravitating 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 G. Voisin et al.: An improved test of the strong equivalence principle with the pulsar in a triple star system different approaches is given in Section 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-scalartensor 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 [START_REF] Capano | GW170817: Stringent constraints on neutron-star radii from multimessenger observations and nuclear theory Carson[END_REF] 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 β 0 -1. In view of this, our limits can be considered as conservative.

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 2018b, 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 equation ( 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 equation ( 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 β 0 , however, gravitational wave tests with binary pulsars become more constraining, in particular for small α 0 . 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 β 0 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. 2019).

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 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). become (almost) independent of the parameter ζ in the sense that ζ s 2 p ∼ 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 equation (35). More generally, in situations where only the strong field of a neutron star can source additional (longrange) 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).

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 [START_REF] Will | Theory and experiment in gravitational physics[END_REF]Damour & Taylor 1992). The mEIH equations of motion describe the first post-Newtonian dynamics of a n-body system which also contains strongly selfgravitating 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 [START_REF] Nordtvedt | [END_REF]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. ). The bodydependent 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 γab , and nine different βa cb . In GR, due to the fulfilment of SEP and the corresponding effacement of the internal structure (see e.g. [START_REF] Damour | The problem of motion in Newtonian and Einsteinian gravity[END_REF]), one has G ab = G N and γab = βa bc = 0.

One can then use the Euler-Lagrange equations [START_REF] Will | Theory and experiment in gravitational physics[END_REF]) to derive the equations of motion for each body: In the weak-field limit one can check that this equation does give the PPN equation of motion (e.g. [START_REF] Soffel | Relativity in astrometry, celestial mechanics, and geodesy[END_REF][START_REF] Will | Theory and experiment in gravitational physics[END_REF].

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 equation (A.1) is derived using the Legendre transform H = a v a • ∂L The momentum of the centre of mass is given by the same expression as in GR only with the replacement G → G ab , P = a ∂L ∂v a and

P = a         m a v a         1 + v 2 a 2c 2 - 1 2 b a G ab m b c 2 r ab         - 1 2 b a G ab m b c 2 r ab (v a • n ab )n ab         .
(A.4)

The centre-of-mass position X satisfies (H /c 2 ) dX dt = P (see e.g. Will 2014b),

Fig. 1 .

 1 Fig. 1. (top) Template pulse profile of PSR J0337+1715 used for timing. 450 hours 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.

Fig. 3 .

 3 Fig. 3. 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. 4 .

 4 Fig. 4. 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.

  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, , 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[START_REF] Ransom | [END_REF]) 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.

Fig. 5 .

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

Fig. 6 .

 6 Fig. 6. 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. Righthand side: distribution of distance to GR derived from the left-hand-side distribution.

Fig. 7 .

 7 Fig. 7. Residuals of the times of arrival using the mean parameters reported in Table2.

Fig. 8 .

 8 Fig. 8. Residuals of the times of arrival corresponding to the parameters of Table2plotted versus the inner-binary orbital phase (top), the outerbinary 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.

Fig. 9 .

 9 Fig. 9. Lomb-Scargle periodogram of the timing residuals for the solution of Table 2. The periodogram is sampled at f m /5 where f m 2268days -1 is the inverse of the full time span. From left to right vertical lines show the frequencies of the red-noise component ( f RN 1650day -1 , black), the Earth orbital period ( f E , green), the outer-orbit orbital period ( f O , orange), the inner-binary orbital period ( f I , red) and the SEP violation signature (2 f If O , purple).

Fig. 10 .

 10 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;Abbott et al. 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 (seeLattimer & Prakash (2001) and https://compose.obspm.fr). Our stiffest EOS, BSk22, is also in agreement with the (more model dependent) constraint of M max 2.3 M byRezzolla 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).

  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 γab = γ ab -1, βa bc = β a bc -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 strongfield parameters satisfy the symmetries G ab = G ba (a b), γab = γba (a b), and βa bc = βa cb (a b, a c

  • n ab ) 2γab (v a -v b ) 2 v b -v a ) n ab • (4v a -3v b -2γ ab (v b -v a )) c 2 n ab 4 + 2γ ac + 2 βa bc .(A.2)

  v a • n ab )(v b • n ab )

Table 1 .

 1 Gaussian priors adopted in MCMC posterior inference.

	Parameter Mean	Std dev	Source
	α	3 h 37 m 43 s .8270	0.37 mas †	1
	δ	17 • 15 14 .8178 0.38 mas †	1
	d	1.3 kpc	160 pc †	3
	µ α

Table 2 .

 2 Voisin et al.: An improved test of the strong equivalence principle with the pulsar in a triple star system Mean values of the MCMC fit with their 68% confidence interval.

	Parameter	Symbol
	Fixed values
	Reference epoch (MJD)	T ref	56492
	Position epoch (MJD)	T pos	57205
	Fitted values
	Right ascension	α	3 h 37 m 43 s .82703(92) +72 -74
	Declination Distance (kpc) Right-ascension proper motion (mas yr -1 ) Declination proper motion (mas yr -1 ) Radial proper motion (mas yr -1 ) Dispersion measure (pc cm -3 ) Dispersion measure variation (pc cm -3 yr -1 )	δ d µ α µ δ µ d DM DM	17 • 15 14 . 818(43) +36 -35 1.3(75) +85 -79 5.(01) +17 -18 -(0.85) +0.69 -0.66 4.(58) +30 -30 21.316(49) +15 -15 (-6) +34 -34 × 10 -6
	Rescaled spin frequency (Hz) Rescaled spin frequency derivative (10 -15 Hz s -1 )	f	f	365.9533379022(28) +24 -25 -2.35587(76) +70 -69
	Orbit of pulsar around centre of mass (CM) of inner binary
	Orbital period (days) Projected semi-major axis (lt-s)	P I a p sin i I	1.62940(06) +30 -30 1.21752(80) +11 -11
	Inclination offset ( • ) Laplace-Lagrange Laplace-Lagrange Time of ascending node (MJD) Long. of asc. nodes offset ( • )	δi = i I -i O e I sin ω I e I cos ω I t ascI δΩ = Ω O -Ω I 0.00001(74) +34 -0.00(29) +29 -29 6.93(65) +19 -19 × 10 -4 -8.5(44) +95 -96 × 10 -5 55917.15(84) +10 -10 -34
	Orbit of CM of inner binary around CM of the whole system
	Orbital period (days) Projected semi-major axis (lt-s) Co-projected semi-major axis (lt-s)	P O a b sin i O a b cos i O	327.255(39) +51 -51 74.6723(74) +57 -58 91.4(35) +47 -48
	Laplace-Lagrange Laplace-Lagrange	e O sin ω O e O cos ω O	0.035114(31) +76 -76 -0.003524(80) +24 -24
	Time of ascending node (MJD)	t ascO	56230.195(11) +41 -41
	Longitude of outer ascending node ( • )	Ω O	-44.(34) +14 -13
	Inner mass ratio SEP ∆ ToA uncertainty rescaling	m i /m p ∆ EFAC	0.1373(50) +18 -18 (4.8) +9.5 -9.4 × 10 -7 1.31(53) +94 -96
	Derived values
	Parallel proper motion (km s -1 ) Plane-of-sky proper motion (km s -1 )	V V ⊥	29.(82) +93 -93 3(3.1) +2.8 -2.5
	Spin frequency (Hz) Spin frequency derivative (10 -15 Hz s -1 )	f f	365.9533630(00) +16 -16 -2.32(44) +41 -37
	Orbit of pulsar around CM of inner binary
	Semi-major axis (lt-s) Orbital inclination ( • )	a p i I	1.924(27) +60 -61 39.2(51) +14 -14
	Orbital eccentricity Longitude of periastron ( • ) Time of periastron passage (MJD) Longitude of asc. node ( • )	e I ω I t p I Ω I	6.98(90) +30 -31 × 10 -4 97.0(22) +76 -75 55917.59(75) +14 -14 -44.(34) +14 -13
	Orbit of CM of inner binary around CM of the whole system
	Semi-major axis (lt-s) Orbital inclination ( • )	a b i O	118.0(53) +36 -37 39.2(37) +14 -14
	Orbital eccentricity Longitude of periastron ( • )	e O ω O	0.035290(78) +78 -78 95.732(19) +27 -27
	Time of periastron passage (MJD) Pulsar mass (M ) Inner-companion mass (M ) Outer-companion mass (M )	t p O m p m i m o	56317.219(76) +52 -52 1.44(01) +15 -15 0.197(80) +19 -19 0.410(58) +40 -40

Article number, page 8 of 22 G.

Nordström's conformally-flat scalar theory, which is also a metric theory, also fulfils the SEP, however, this is excluded by Solar System experiments (Deruelle

2011).

Shao (2016) investigates the possibility of constraining a difference in active and passive gravitational mass with the pulsar system under consideration in this paper.

More precisely, the constraint is on ∆ pulsar,Galaxy -∆ companion,Galaxy .

tempo is a standard pulsar timing software, this can be found at http://tempo.sourceforge.net.

http://psrchive.sourceforge.net.

Preferred location effects are already tightly constrained using pulsars byShao & Wex (2013).

7 Source, data, and results are available at: doi.org/10.5281/ zenodo.3778978

Although addition of velocities only applies to Newtonian mechanics, we are here only interested in transforming the orbital elements into initial conditions for the numerical integrator. Such transformation is somewhat arbitrary, and we choose to add these velocities for simplic-

ity.9 Boost library version 1.55.0 www.boost.org

Article number, page 10 of 22 G.Voisin et al.: An improved test of the strong equivalence principle with the pulsar in a triple star system

The number between brackets gives the uncertainty on the last digit(s).

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

Article number, page 20 of 22 G. Voisin et al.: An improved test of the strong equivalence principle with the pulsar in a triple star system 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.