© The Authors 2024

1 Introduction

Planetary ephemerides have evolved over the years in line with improvements in the observational accuracy of the astrometry of planets and natural satellites using the navigation tracking of these systems. The motion of the celestial bodies in our Solar System can be solved directly by integrating their equations of motion numerically. The improved accuracy in measurements of observables from space missions (e.g. Cassini-Huygens, Mars Express, BepiColombo, Juice, etc.) makes the Solar System a suitable environment for testing general relativity theory (GRT) as well as alternative theories of gravity by means of Solar System ephemerides. For a detailed review regarding the use of planetary ephemerides for testing gravity, we refer to Fienga & Minazzoli (2024). Since 2003, the INPOP (Intégrateur Numérique Planétaire de l’Observatoire de Paris) planetary ephemerides have been developed (Fienga et al. 2008, 2009) with numerical integration of the Einstein-Infeld-Hoffmann equations of motion proposed by Moyer (2003), and by fitting the parameters of the dynamical model to the most accurate planetary observations following Soffel et al. (2003). In the present work, we use the INPOP21a planetary ephemerides (Fienga et al. 2021), which benefit from the latest Juno and Mars orbiter tracking data up to 2020, and from radio-science observations from space missions such as Cassini-Huygens, Mars Express, and Venus Express. Moreover, a fit of the Moon–Earth system to Lunar Laser Ranging (LLR) observations up to 2020 is included. The dynamical modelling used in INPOP21a takes into account the eight planets of the Solar System, as well as Pluto, the Moon, the asteroids, the Kuiper belt, the Sun oblateness, and the Lense–Thirring effect. We refer the reader to Fienga et al. (2021) for a detailed review of this version of INPOP. In recent years, the post-Newtonian (PPN) parameters have already been tested with planetary ephemerides (Fienga & Minazzoli 2024). However, previous derivations did not take into account the potential violation of the strong equivalence principle (SEP) owing to its assumed small impact on the ephemerides. In other words, the effect of the gravitational self-energy of the planets violating the SEP was not implemented in the equations of motion. This decision was based on order-of-magnitude derivations that indicated that the effect was not deemed significant given the accuracy of the previous ephemerides.

Nevertheless, almost all alternative theories predict a violation of the SEP – that is, they predict that bodies with different self-gravitating energies fall differently in a given gravitational field (Nordtvedt 1968a,b; Will 2018). In the PPN framework, this effect is parameterized by the Nordtvedt parameter η – see Eq. (6) – in reference to Nordtvedt (1968a,b), and η = 0 in general relativity. Therefore, it was only a matter of time before planetary ephemerides reached a level of accuracy that necessitated that this effect be accounted for in the integration of the equations of motion.

In the present work, we implement the whole set of relevant equations in the Brans–Dicke case – including the effects of the predicted violation of the SEP. We use the recent work of Bernus et al. (2022), who derived the whole phenomenology of Einstein-dilaton theories, and restrict our analysis to the one-parameter family that corresponds to Brans–Dicke theories. As we show in Sect. 2.1, one of the reasons for caution is that there is a potential ambiguity regarding the types of mass (gravitational or inertial) that appear in various equations, namely the equation of motion, the Shapiro delay, and the definition of the barycenter. We therefore start from a consistently derived phenomenology.

The motivation for investigating violations of the SEP in planetary ephemerides within the relatively limited framework of Brans–Dicke theories is two-fold. Firstly, as explained in a recent review (Fienga & Minazzoli 2024), because of the correlations between the parameters of the Solar System (masses, initial conditions, shape, etc.) and the theoretical parameters that parameterize deviations from the laws of general relativity (e.g. PPN parameters), it is somewhat difficult to produce statistically robust tests of alternative theories to general relativity. Fienga & Minazzoli (2024) indeed notably argue that one should focus on each set of theories in a consistent manner, and this is our aim for the present study. Secondly, a robust statistical Bayesian methodology was recently developed and used to test a potential Yukawa suppression of the Newtonian potential in the Solar System (Mariani et al. 2023b). It is interesting to apply this new methodology to other alternative phenomenologies. However, it is only feasible to apply such an approach to phenomenologies that have just one theoretical parameter to be tested. Therefore, if the methodology is applied to scalar-tensor theories, one must restrict the study to the Brans–Dicke case.

2 Dynamical and statistical formalism

2.1 Gravitational framework

A discrepancy in the free fall of particles with differing compositions in a particular gravitational field would signify a violation of the weak equivalence principle (WEP). Meanwhile, a deviation for extended bodies with distinct gravitational self-energies would indicate a breach of the gravitational weak equivalence principle (GWEP). The latter forms a component of the SEP, as delineated by Will (2018). In all instances, through planetary ephemerides, our objective is to verify the universality of free fall (UFF), or, in other words, to confirm whether all bodies, regardless of their characteristics, experience identical rates of fall in a given gravitational field.

The violations of both the WEP and GWEP share a characteristic feature: at the Newtonian level, their equations of motion are expressed as follows (Bernus et al. 2022; Minazzoli & Hees 2016): $$a_{\text{T}}=-{\displaystyle \sum _{\text{A}\ne \text{T}}{\mu }_{\text{A}}}\frac{r_{\text{AT}}}{r_{\text{AT}}^3}\left(1+{\delta }_{\text{T}}+{\delta }_{\text{AT}}\right),$$(1)

where a_T is the coordinate acceleration of body T (whereas the letter A is used as index for the other bodies); r_AT is the relative position of body T with respect to A; r_AT = |r_AT|, with ${\mu }_{\text{A}}=G\times m_{\text{A}}^G$ being the body A gravitational parameter; and δ_T and δ_AT are coefficients that depend on the composition of the bodies T and A for the case of the violation of the WEP only. The most general post-Newtonian phenomenology for Einstein-dilaton theories was derived by Bernus et al. (2022), and depends on several parameters that are related to the coupling between matter and gravity. Because the exploration of a high-dimensional theoretical space is very demanding in terms of computation time, a reduction of the number of parameters was made owing to the roughly similar compositions of the celestial bodies, and a simple rejection sampling was applied to the solutions of the ephemerides in order to estimate the constraints of the theoretical parameters.

In what follows, we use an opposing strategy: given the recent very effective Bayesian methodology that is well suited to the one-parameter family of theories, we reduce the general Einstein-dilaton action to the case of only one parameter, effectively recovering the Brans–Dicke theory, and we apply the Bayesian methodology to this case.

Following Bernus et al. (2022), the equation of motion of the Brans–Dicke special case corresponds to $$\begin{array}{l}{a}_{\text{T}}=-{\displaystyle \sum _{\text{A}\ne \text{T}}\frac{{\mu }_{\text{A}}}{r_{\text{AT}}^3}}{r}_{\text{AT}}\left(1+{\delta }_{\text{T}}\right)-{\displaystyle \sum _{\text{A}\ne \text{T}}\frac{{\mu }_{\text{A}}}{r_{\text{AT}}^3{c}^2}}{r}_{\text{AT}}\{\gamma v_{\text{T}}^2+(\gamma +1){\upsilon }_{\text{A}}^2\\ -2(1+\gamma ){\upsilon }_{\text{A}}\cdot {\upsilon }_{\text{T}}-\frac{3}{2}{\left(\frac{r_{\text{AT}}\cdot {\upsilon }_{\text{A}}}{r_{\text{AT}}}\right)}^2-\frac{1}{2}{r}_{\text{AT}}\cdot a_{\text{A}}-2\gamma {\displaystyle \sum _{B\ne T}\frac{{\mu }_{\text{B}}}{r_{TB}}}+{\displaystyle \sum _{B\ne A}\frac{{\mu }_{\text{B}}}{r_{\text{AB}}}}\}\\ +{\displaystyle \sum _{\text{A}\ne \text{T}}\frac{{\mu }_{\text{A}}}{c^2{r}_{\text{AT}}^3}}\left[2(1+\gamma )r_{\text{AT}}\cdot {\upsilon }_{\text{T}}-(1+2\gamma )r_{\text{AT}}\cdot {\upsilon }_{\text{A}}\right]\\ \left({\upsilon }_{\text{T}}-{\upsilon }_{\text{A}}\right)+\frac{3+4\gamma }{2}{\displaystyle \sum _{\text{A}\ne \text{T}}\frac{{\mu }_{\text{A}}}{c^2{r}_{\text{AT}}}}{a}_{\text{A}},\end{array}$$(2)

where δ_T and the parameter γ depend only on a universal coupling constant α₀, such that $$\gamma =\frac{\left(1-{\alpha }_0^2\right)}{\left(1+{\alpha }_0^2\right)}$$(3)

and $${\delta }_{\text{T}}=-(1-\gamma ){\text{Ω}}_{\text{T}}\hspace{0.17em}:=\eta {\text{Ω}}_{\text{T}},$$(4)

with $${\text{Ω}}_{\text{T}}=\frac{3}{5}\frac{{\mu }_{\text{T}}}{c^2{R}_{\text{T}}},$$(5)

where R_T is the radius of the planet T. This implies that the Nordtvedt parameter in this case is $$\eta =-(1-\gamma )=-2\frac{{\alpha }_0^2}{1+{\alpha }_0^2}.$$(6)

The Shapiro delay is also modified and is now such that the whole coordinate propagation time between an emitter and a receiver (t_r − t_e) reads as follows: $$c\left(t_r-t_e\right)=R+{\displaystyle \sum _{\text{T}}\left(1+\gamma -{\delta }_{\text{T}}\right)}\frac{{\mu }_{\text{T}}}{c^2}\ln \frac{n\hspace{0.17em}\cdot \hspace{0.17em}{r}_{rT}\hspace{0.17em}+\hspace{0.17em}{r}_{rT}\hspace{0.17em}+\hspace{0.17em}\frac{4{\mu }_{\text{T}}}{c^2}}{n\hspace{0.17em}\cdot \hspace{0.17em}{r}_{eT}\hspace{0.17em}+\hspace{0.17em}{r}_{eT}\hspace{0.17em}+\hspace{0.17em}\frac{4{\mu }_{\text{T}}}{c^2}},$$(7)

where R is the coordinate Euclidean distance between the emitter and the receiver. Another way to formulate the Shapiro delay is in terms of the inertial gravitational parameter, $${\mu }_{\text{T}}^i=G\times m_{\text{T}}^i=\left(1-{\delta }_{\text{T}}\right){\mu }_{\text{T}},$$(8)

where $m_{\text{T}}^i$ is the inertial mass of the body T, such that $$c\left(t_r-t_e\right)=R+{\displaystyle \sum _{\text{T}}(1+\gamma )}\frac{{\mu }_{\text{T}}^i}{c^2}\ln \frac{n\hspace{0.17em}\cdot \hspace{0.17em}{r}_{rT}\hspace{0.17em}+\hspace{0.17em}{r}_{rT}\hspace{0.17em}+\hspace{0.17em}\frac{4{\mu }_{\text{T}}}{c^2}}{n\hspace{0.17em}\cdot \hspace{0.17em}{r}_{eT}\hspace{0.17em}+\hspace{0.17em}{r}_{eT}\hspace{0.17em}+\hspace{0.17em}\frac{4{\mu }_{\text{T}}}{c^2}}.$$(9)

In other words, the Shapiro delay is sensitive to the iner-tial masses, whereas the equation of motion depends on the gravitational masses.

Equally importantly, the definition of the Solar System barycenter (SSB) also involves the inertial gravitational parameter ${\mu }_{\text{A}}^i$. Indeed, from the Lagrangian formulation of the equations of motion, Bernus et al. (2022) show that the following barycenter constant vector is a first integral of the equations of motion, $$q=G-Vt\text{,}$$(10)

where $$G=\frac{c^2}{h}{\displaystyle \sum _{\text{A}}{\mu }_{\text{A}}^i}{z}_{\text{A}}\left(1+\frac{v_{\text{A}}^2}{2c^2}-\frac{1}{2c^2}{\displaystyle \sum _{B\ne A}\frac{{\mu }_{\text{B}}}{r_{\text{AB}}}}\right)$$(11)

are the coordinates of the relativistic barycenter of the system and $$V=\frac{c^{2}P}{h}$$(12)

is the velocity of the barycenter motion. h is the conserved energy – whose value does not affect what follows but can be found in Bernus et al. (2022) – and P is the conserved linear momentum, which reads $$\begin{array}{ll}\hfill & \underset{\text{A}}{{\displaystyle P=\sum }}{p}_{\text{A}}\hfill \\ \hfill & \underset{\text{A}}{{\displaystyle =\sum }}{\mu }_{\text{A}}^i{\upsilon }_{\text{A}}\left[1+\frac{1}{2c^2}\left({\upsilon }_{\text{A}}^2-{\displaystyle \sum }_{B\ne A}\frac{{\mu }_{\text{B}}}{r_{\text{AB}}}\right)\right]\hfill \\ \hfill & -\frac{1}{2c^2}{\displaystyle \sum }_{\text{A}}{\displaystyle \sum }_{B\ne A}\frac{{\mu }_{\text{A}}{\mu }_{\text{B}}}{r_{\text{AB}}}\left(n_{\text{AB}}\cdot {\upsilon }_{\text{A}}\right)n_{\text{AB}}.\hfill \end{array}$$(13)

In this work, we implemented Eqs. (2) and (7) in the INPOP planetary ephemerides numerical integration of the planet equations of motion and its adjustment to observations, whereas – because the system has a conserved linear momentum – the equation for the SSB (defined with the inertial gravitational parameter) is integrated once at J2000, which is the initial date of integration for INPOP ephemerides. We used INPOP21a for this test, as in Mariani et al. (2023b,a), by considering different values of α₀ and the changes this induces in γ and η. The obtained results are presented in Sect. 3.2 in terms of γ, while Sect. 2.2 explains the statistical method.

To assess the effect of considering the potential violation of the GWEP (and therefore of the SEP) in planetary ephemerides, we constrained the value of γ either assuming that η = 0 (in Sect. 3.1) – as was the case in previous planetary ephemerides – or assuming that η = −(1 − γ) (in Sect. 3.2).

2.2 Statistical aspects

Starting from Fienga et al. (2015; also in Di Ruscio et al. 2020 and Fienga et al. 2020), several tools have been used to determine the goodness of the INPOP fit with respect to modifications in the equations of motion or in the framework of the ephemeris, as well as the observations. Within this context, the computation of the INPOP χ² plays a role in determining which model (or which data) leads to a significant improvement of the ephemeris computation. In the present work, the computation of χ² is the output – for a given value of γ – of the INPOP iterative fit, that is, after the adjustment of all its astronomical parameters (see Mariani et al. 2023b,a; Bernus et al. 2019). The χ²(γ) is computed following (14) $${\chi }^2(\gamma ,k)\equiv \frac{1}{N_{\text{obs}}}{\displaystyle \sum }_{i=1}^{N_{\text{obs}}}{\left(\frac{g^i(\gamma ,k)-d_{\text{obs}}^i}{{\sigma }_i}\right)}^2,$$(14)

where γ is a fixed value, k are the astronomical parameters fitted with INPOP, N_obs is the number of observations, the function ɡⁱ represents the computation of observables, the vector $d_{\text{obs}}={\left(d_{\text{obs}}^i\right)}_i$ is the vector of observations and σ_i are the observational uncertainties. The goal of the methodology applied is to obtain a posterior for the parameter γ. The general pipeline used to obtain such a posterior is similar to that adopted by Mariani et al. (2023a,b). We compute the value χ²(γ) (see Eq. (14)) for several different values of γ, spreading over our domain of interest. For any given γ, the value χ²(γ) is obtained as an outcome of the full INPOP iterative adjustment, with γ fixed. In this way, the astronomical parameters are adjusted with the least-squares procedure, whereas γ is a fixed value. Fixing the value of γ is necessary as it ensures the γ contribution to the dynamics, avoiding high correlations between γ and the other astronomical parameters k (see Anderson et al. 1978; Bernus et al. 2019) for further details about correlation problem). In order to produce the posterior, we use a Markov chain Monte Carlo (MCMC) method, the Metropolis–Hastings (MH) algorithm. In order to produce the MCMC, it is necessary to evaluate the likelihood, sequentially, and therefore the χ², for thousands of different values of γ (see Mariani et al. 2023a for further details). As the computation of χ²(γ) is expensive in terms of time, it becomes difficult to run the MCMC with the χ² direct calculation. In order to overcome the problem of computation time, we apply a Gaussian process regression (GPR) to interpolate among the values (γ, χ²(γ)) already computed beforehand. Starting from this set of points, we obtain, thanks to the GPR, a function $\gamma \mapsto {\tilde{\chi }}^2(\gamma )$, together with an uncertainty relative to the possible error of interpolation. The regression $\gamma \mapsto {\tilde{\chi }}^2(\gamma )$ is necessary to run the MH algorithm, which has, as an outcome, the posterior density. Moreover, an a priori density distribution has to be provided for the parameter of interest (i.e. γ) to sample with the MCMC from the posterior. We used a uniform prior distribution for γ. The MH algorithm is based on the idea that the drawings from the posterior distribution are done sequentially and according to a random acceptance process. The outcome of the algorithm is a sequence of random samples from the posterior (this sequence is the Markov chain) and its equilibrium distribution is the posterior itself. At each step of the algorithm, one new element of the sequence is computed and proposed, and is accepted or rejected. The acceptance/rejection step is random: this randomness is determined by a certain probability (which changes at each step) based on the last accepted element of the sequence as well as the likelihood of the last element and of the new candidate element.

3 Results

3.1 Limit on γ without SEP violation

In order to evaluate the impact of the method within a well-known context Mariani et al. (2023a), we first estimated the possible violation of general relativity by considering 1 − γ ≠ 1 without considering possible SEP violation (i.e. η = 0). In this case, the terms δ_T in Eq. (2) and δ_A in Eq. (7) vanish and only the terms with γ remain. With this configuration, our results are comparable to classical conjunction tests like the one obtained during the Cassini interplanetary phase when the radio-science signal of the Cassini s/c grazed the Sun before reception on earth. During this conjunction, the Shapiro delay is at its maximum and Bertotti et al. (2003) deduced the best constraint on γ so far with a global fit of s/c orbit parameters and γ. As this part of the mission (interplanetary phase) is supposed to be the least affected by complex gravitational interactions or maneuvers, the determination is in principle the most accurate and the least affected by biases or noise. Because of the estimation method (least square), the Bertotti et al. (2003) constraint, 1 − γ = (−2.1 ± 2.3) × 10⁻⁵, is given at 1σ, assuming a Gaussian distribution of the noise. Two main differences between our approach and the results obtained by Bertotti et al. (2003) have to be stressed. First, by construction of our model, only values of γ where 1 − γ > 0 are tested. Once we compare our results with Bertotti et al. (2003), we have to compare with the absolute value of their interval as it is represented in Fig. 1. Second, we use a Bayesian approach using a uniform prior with a maximum value for 1 − γ of 15 × 10⁻⁵ and providing a posterior distribution of acceptable 1 − γ values, as discussed in Sect. 2.2. In this context, for comparisons of results, it is important to consider the confidence level (C.L.) deduced from the posteriors. In Fig. 1 and Table 1, we see the limit of about 4.4 × 10⁻⁵ at 66.7% C.L. (1σ) of the ‖1 − γ‖ interval deduced from Bertotti et al. (2003; black dashed line on Fig. 1) and in yellow, the posterior distribution of ‖1 − γ‖ obtained with this work. At 99.7% C.L., we obtain a limit of 2.5 × 10⁻⁵ (red dot-dashed line on Fig. 1) and of about 1.73 × 10⁻⁵ at 66.7% C.L. (blue dotted line on Fig. 1). It is interesting to note, at this stage, that the limits we obtain are of the same order of magnitude as those obtained by Bertotti et al. (2003). However, we stress that these results were obtained while exploring only γ values smaller than 1.

Fig. 1 Posterior distributions without SEP violation. Yellow: Posterior distribution for 1 − γ without the implementation of the SEP violation. Green: Absolute value of the simulated posterior for the determination of 1 − γ with Cassini. Red line (dot-dashed): 99.7% C.L. ((1 − γ) = 2.50 × 10−5). Blue line (dotted): 66.7% C.L. ((1 − γ) = 1.73 × 10−5). Purple line (dashed): 95% C.L. ((1 − γ) = 2.40 × 10−5). The prior for (1 − γ) is a uniform prior between 0 and 15 × 10−5. Solid black line: Absolute value of Cassini central value estimation ((1 − γ) = −2.1 × 10−5). Dashed black line: Upper bound estimation from Cassini at 1σ C.L.

3.2 Limit on γ and α₀ with SEP violation

After considering the results of an exploration of the 1 − γ without SEP violation, we now switch to the full Brans–Dicke formalism with SEP violation as described in Sect. 2.1. We used the method described in Sect. 2.2 for the determination of acceptable values for the parameter α₀. Figure 2 shows the posterior of 1 − γ and the deduced α₀ from Eq. (3). At 99.7 % C.L., we obtain a limit on ||1 − γ|| of about 2.83 × 10⁻⁵, leading to a constraint on α₀ of about 3.76 × 10⁻³. From the constraint on α₀ and on γ, one can deduce limits for the Nordtvedt parameter η = γ − 1. As given on Table 1, at 99.7 % C.L., η > −2.83 × 10⁻⁵, at 95 % C.L., η > −2.70 × 10⁻⁵, and at 66.7 % C.L., η > −1.92 × 10⁻⁵. We stress that these constraints are only for η < 0 (see Eq. (6)). Different values of η can be compared with these results (see Table 1) but all are deduced from a direct fit either of the s/c orbit, such as with Genova et al. (2018), or of the Earth–Moon system (Viswanathan et al. 2018).

In Genova et al. (2018), η is obtained during the orbit determination of the s/c Messenger orbiting Mercury. In the global fit, in addition to orbit determination parameters, the PPN parameter β, the Sun oblateness J₂, and η are also estimated. PPN γ is fixed to the Bertotti et al. (2003) value.

As discussed in Fienga & Minazzoli (2024), the four parameters are strongly correlated, with or without the use of the Nordtvedt relation, $$\eta =4(\beta -1)-(\gamma -1).$$(15)

For comparison, in Table 1, we propose some values found in the literature obtained with very different methods. In particular, we consider the thresholds obtained with Bernus et al. (2022) in the context of the Einstein-dilaton theories, which can be seen as a generalisation of the Brans–Dicke scalar–tensor theories. Indeed, in this case, coupling between matter and the scalar field is introduced and can take either a linear or a non-linear form. In Bernus et al. (2022), limits on the universal coupling α₀ have been obtained using planetary ephemerides (INPOP19a) together with constraints for the gaseous and the telluric planets, α_G and α_T, respectively. The comparison between the limit obtained for α₀ in this work and that from Bernus et al. (2022) is not straight-forward, as in our case, we have α_G = α_T = 0.

4 Discussion

4.1 The equivalence principle

It is interesting to note that the limit obtained for 1 − γ in the BD framework with the SEP violation is larger than that obtained without the SEP violation (see Fig. 1). We interpret this as marginal evidence that, at the level of accuracy of current planetary ephemerides, the effect of the SEP violation is no longer negligible. Indeed, adding the effect increases the correlations between γ and planetary parameters, which, here, comes from the additional correlations that exist between δ_T and planetary parameters (masses, initial conditions, etc.). This leads to the usual reduction in the accuracy of the estimation of the constraint on the theoretical parameter being tested. Indeed, with a greater number of correlations, more freedom is provided with which to better fit the data; this leads to weaker constraints. For a comprehensive discussion about the general issue of correlations between theoretical and planetary parameters, we refer the reader to Fienga & Minazzoli (2024). Also, because we observe that including the SEP violation leads to slightly more stable solutions for the planetary orbits, we believe that including the SEP violation in alternative theories – that is, η ≠ 0 – is more consistent at the dynamical level than not including it.

Fig. 2 Posterior distribution for 1 − γ and the corresponding values of α0 in the Brans–Dicke framework with the implementation of the SEP violation. Red line (dot-dashed): 99.7% C.L. ((1 − γ) = 2.83 × 10−5 and α0 = 3.762 × 10−3). Blue line (dotted): 66.7% C.L. ((1 − γ) = 1.92 × 10−5 and α0 = 3.096 × 10−3). Purple line (dashed): 95% C.L. ((1 − γ) = 2.70 × 10−5 and α0 = 3.675 × 10−3). The prior for (1 − γ) is a uniform prior between 0 and 15 × 10−5.

4.2 Conversion in terms of the usual Brans–Dicke parameter ω_BD

One can simply convert the constraint on α₀ to a constraint on the usual Brans–Dicke parameter ω_BD (Will 2014) $${\omega }_{\text{BD}}=\frac{1}{2{\alpha }_0^2}-\frac{3}{2}.$$(16)

Therefore, based on the results given in Table 1, we can derive a 99.7% C.L. with ω_BD > 35 382. The limits obtained with the planetary ephemerides in this work are not as stringent as those obtained by Voisin et al. (2020). Indeed, using observations of a pulsar in a triple-star system, Voisin et al. (2020) deduce that ω_BD > 130000 although this value depends on the equation of state being assumed for the neutron star, and variation of a few tens of percent is generally found when using different equations of state (Voisin et al. 2020). Hence, constraints from planetary ephemerides are currently less powerful than the best constraints from pulsars. Nevertheless, we anticipate that future spacecraft missions, such as BepiColombo, will significantly improve the constraining power of planetary ephemerides in the near future.

5 Conclusions

In the present paper, we used a recent powerful Bayesian methodolgy developed for a one-parameter theoretical model in order to obtain the latest constraints from planetary ephemerides on the Brans–Dicke theory of gravity. In this context, we find that ‖1 − γ‖ < 1.92 ×10⁻⁵ at the 66.7% C.L., while the previous best Solar System constraints from ranging data of the Cassini spacecraft (Bertotti et al. 2003) on γ led to ‖1 − γ‖ < 4.4 × 10⁻⁵ at 66.7% C.L. At the 99.7% C.L., the new planetary ephemerides constraint reads ‖1 − γ‖ < 2.83 × 10⁻⁵. While this is still inferior to the best constraints from pulsars (Voisin et al. 2020), missions such as BepiColombo should improve this constraint in the near future. We also report marginal evidence suggesting that the effect of the violation of the strong equivalence principle can no longer be assumed to be negligible with the current accuracy of planetary ephemerides.

Acknowledgements

V.M. was funded by CNES (French Space Agency) and UCA EUR Spectrum doctoral fellowship. This work was supported by the French government, through the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR) under reference number ANR-15-IDEX-01. The authors are grateful to the OPAL infrastructure and the Université Côte d’Azur Center for High-Performance Computing for providing resources and support. The authors thank G. Voisin for his useful inputs and discussions. This work was supported by the French group for gravitation and metrology PNGRAM (Programme National Gravitation, Références, Astronomie, Métrologie).

References

All Tables

All Figures

Fig. 1 Posterior distributions without SEP violation. Yellow: Posterior distribution for 1 − γ without the implementation of the SEP violation. Green: Absolute value of the simulated posterior for the determination of 1 − γ with Cassini. Red line (dot-dashed): 99.7% C.L. ((1 − γ) = 2.50 × 10−5). Blue line (dotted): 66.7% C.L. ((1 − γ) = 1.73 × 10−5). Purple line (dashed): 95% C.L. ((1 − γ) = 2.40 × 10−5). The prior for (1 − γ) is a uniform prior between 0 and 15 × 10−5. Solid black line: Absolute value of Cassini central value estimation ((1 − γ) = −2.1 × 10−5). Dashed black line: Upper bound estimation from Cassini at 1σ C.L.

In the text

Fig. 2 Posterior distribution for 1 − γ and the corresponding values of α0 in the Brans–Dicke framework with the implementation of the SEP violation. Red line (dot-dashed): 99.7% C.L. ((1 − γ) = 2.83 × 10−5 and α0 = 3.762 × 10−3). Blue line (dotted): 66.7% C.L. ((1 − γ) = 1.92 × 10−5 and α0 = 3.096 × 10−3). Purple line (dashed): 95% C.L. ((1 − γ) = 2.70 × 10−5 and α0 = 3.675 × 10−3). The prior for (1 − γ) is a uniform prior between 0 and 15 × 10−5.

In the text