Likewise, at the level of the solar system barycenter, the TDB is defined as a linear transformation of the TCB. The relations between the various time scales is shown in Fig. 5.

The difference between TT and TDB is produced by planetary ephemerides, by integrating the following equation together with the equations of motion (Klioner 2008; Fienga et al. 2009)

$$ \frac{d(TT - TDB)}{d(TDB)} = \left(L_B + \frac{1}{c_2}a\right)(1 + L_B - L_G) - L_G + \frac{1}{c_4}b, \quad (85) $$

with $L_B$ and $L_G$ are defining constants for TDB relatively to TCB and TT relatively to TCG, respectively (see e.g. Klioner 2008, Petit and Luzum 2010 for the full definition) and where

$$ \begin{aligned} a &= -\frac{1}{2}v_T^2 - \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \\ b &= -\frac{1}{8}v_T^4 + \frac{1}{2} \left[ \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \right]^2 + \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \left\{ 4v_T \cdot v_A - \frac{3}{2}v_T^2 - 2v_A^2 \right. \\ & \qquad \left. + \frac{1}{2}a_A \cdot r_{AT} + \frac{1}{2} \left( \frac{v_A \cdot r_{AT}}{r_{AT}} \right)^2 + \sum_{B \neq A} \frac{\mu_B}{r_{BA}} \right\}, \end{aligned} \quad (86) $$

following the notations of Eq. (84). Park et al. (2021), terms related to the perturbations induced by the oblateness of the Sun were also added in the DE440 (TT-TDB) computation.

Finally, planetary observations are also related to time. Following the IERS conventions (Petit and Luzum 2010), those observations obtained on the ground are given in UTC, related to TDB by TT and TAI. The observations obtained in other planetary systems (e.g. positions of a spacecraft orbiting Mercury) are, up to now, given also in UTC as the differences between the coordinate time defined for the corresponding planetary system (in the previous example, Mercury) and TDB or TT, are not significant at the present day accuracy. However, with missions such as Bepi-Colombo, it will be necessary to account for the local gravitational potential in the definition of the observational time scales (e.g. Turyshev et al. 2013a, Nelson and Ely 2006, Milani et al. 2002). In conclusion, it is important to keep in mind that the various measured proper times are converted beforehand into coordinate times— according to the laws of general relativity.

### 3.2.3 Inertial frame

The coordinate frame of the planetary and lunar ephemerides is linked to the International Celestial Reference System (ICRS) by its current realization achieved by VLBI measurements of the positions of extragalactic radio sources (e.g., quasars) defined in the ICRF (Standish 1998b). The planetary orbits are tied to ICRF because the observations used for their adjustment have been obtained in the ICRF. For the inner planets, the VLBI observations of Venus and Mars-orbiting missions give a tie
Likewise, at the level of the solar system barycenter, the TDB is defined as a linear transformation of the TCB. The relations between the various time scales is shown in Fig. 5.

The difference between TT and TDB is produced by planetary ephemerides, by integrating the following equation together with the equations of motion (Klioner 2008; Fienga et al. 2009)

$$ \frac{d(TT - TDB)}{d(TDB)} = \left(L_B + \frac{1}{c_2}a\right)(1 + L_B - L_G) - L_G + \frac{1}{c_4}b, \quad (85) $$

with $L_B$ and $L_G$ are defining constants for TDB relatively to TCB and TT relatively to TCG, respectively (see e.g. Klioner 2008, Petit and Luzum 2010 for the full definition) and where

$$ \begin{aligned} a &= -\frac{1}{2}v_T^2 - \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \\ b &= -\frac{1}{8}v_T^4 + \frac{1}{2} \left[ \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \right]^2 + \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \left\{ 4v_T \cdot v_A - \frac{3}{2}v_T^2 - 2v_A^2 \right. \\ & \qquad \left. + \frac{1}{2}a_A \cdot r_{AT} + \frac{1}{2} \left( \frac{v_A \cdot r_{AT}}{r_{AT}} \right)^2 + \sum_{B \neq A} \frac{\mu_B}{r_{BA}} \right\}, \end{aligned} \quad (86) $$

following the notations of Eq. (84). Park et al. (2021), terms related to the perturbations induced by the oblateness of the Sun were also added in the DE440 (TT-TDB) computation.

Finally, planetary observations are also related to time. Following the IERS conventions (Petit and Luzum 2010), those observations obtained on the ground are given in UTC, related to TDB by TT and TAI. The observations obtained in other planetary systems (e.g. positions of a spacecraft orbiting Mercury) are, up to now, given also in UTC as the differences between the coordinate time defined for the corresponding planetary system (in the previous example, Mercury) and TDB or TT, are not significant at the present day accuracy. However, with missions such as Bepi-Colombo, it will be necessary to account for the local gravitational potential in the definition of the observational time scales (e.g. Turyshev et al. 2013a, Nelson and Ely 2006, Milani et al. 2002). In conclusion, it is important to keep in mind that the various measured proper times are converted beforehand into coordinate times— according to the laws of general relativity.

### 3.2.3 Inertial frame

The coordinate frame of the planetary and lunar ephemerides is linked to the International Celestial Reference System (ICRS) by its current realization achieved by VLBI measurements of the positions of extragalactic radio sources (e.g., quasars) defined in the ICRF (Standish 1998b). The planetary orbits are tied to ICRF because the observations used for their adjustment have been obtained in the ICRF. For the inner planets, the VLBI observations of Venus and Mars-orbiting missions give a tie
with an accuracy better than the milliarcsecond (Fienga et al. 2011a; Folkner et al. 2014). For the outer planets, the link is maintained with the same level of uncertainty thanks the VLBI observations of Jones et al. (2019), Fienga et al. (2011a), and Folkner et al. (2014) orbiting Jupiter and Saturn, respectively. Some other methods have been investigated for enhancing the link between planetary ephemerides and ICRF. In particular, one can notice the use of Lunar Laser Ranging observations (Pavlov 2020) and the use of GAIA DR2 asteroid positions (Deram et al. 2022). The tie between the planetary ephemerides frame and ICRF is confirmed to be sub-mas level accuracy in both approaches.

### 3.2.4 Definition of the solar system barycenter (SSB)

The definition of the solar system barycenter at the origin of the time of integration is based on the hypothesis of the conservation of the energy of the dynamical system composed by the planets, the Moon and the asteroids (including Trans-Neptunian Objects) in the solar system. The position $R$ and the velocity $V$ of the center of mass are then invariant. Using the notation of DE430 (Folkner et al. 2014), it comes

$$ \left\{ \begin{aligned} R &= \frac{\sum_i \mu_i^* \mathbf{r}_i}{\sum_i \mu_i^*} = 0 \\ V &= \frac{dR}{dt} = 0 \end{aligned} \right. \qquad (87) $$

where $\mathbf{r}_i$ and $\mathbf{v}_i$ being respectively its barycentric position and velocity vectors of the planet $i$. $\mu^*$ is given by

$$ \left\{ \begin{aligned} \mu_i^* &= \mu_i \left( 1 + \frac{\mathbf{v}_i^2}{2c^2} - \frac{1}{2c^2} \sum_{j \neq i} \frac{\mu_j}{r_{ij}} \right) \\ \dot{\mu}_i^* &= \frac{\mu_i}{2c^2} \left( \sum_{j \neq i} \mu_j \frac{(\mathbf{r}_j - \mathbf{r}_i) \cdot (\mathbf{v}_j + \mathbf{v}_i)}{r_{ij}^3} \right) \end{aligned} \right. \qquad (88) $$

with $\mu_i$ the product of the gravitational constant $G$ with the inertial mass of the body $i$. The term $\dot{\mu}_i^*$ is not included in DE430 nor DE440 ephemerides (Folkner et al. 2014; Park et al. 2021) but is accounted for in the INPOP ephemerides (Fienga et al. 2008). The positions and velocities of the SSB are then obtained by integration of the following equations and up to the $c^{-2}$ order,

$$ \left\{ \begin{aligned} \sum_i \mu_i^* \mathbf{r}_i &= 0 \\ \sum_i \mu_i^* \mathbf{v}_i + \dot{\mu}_i^* \mathbf{r}_i &= 0 \end{aligned} \right. \qquad (89) $$

In the INPOP formalism (Fienga et al. 2008), Eqs. (89) and (88) are solved only at the initial step of the planetary integration and a constant vector in positions and velocities is subtracted to all the body positions and velocities for having $R$ and $V$ to 0. Once the frame is centered on the SSB defined by Eq. (89) at J2000 (initial date of integration for INPOP ephemerides), the equations of motion of the solar system bodies and of the Sun are integrated in this fixed frame. The method described here is
with an accuracy better than the milliarcsecond (Fienga et al. 2011a; Folkner et al. 2014). For the outer planets, the link is maintained with the same level of uncertainty thanks the VLBI observations of Jones et al. (2019), Fienga et al. (2011a), and Folkner et al. (2014) orbiting Jupiter and Saturn, respectively. Some other methods have been investigated for enhancing the link between planetary ephemerides and ICRF. In particular, one can notice the use of Lunar Laser Ranging observations (Pavlov 2020) and the use of GAIA DR2 asteroid positions (Deram et al. 2022). The tie between the planetary ephemerides frame and ICRF is confirmed to be sub-mas level accuracy in both approaches.

### 3.2.4 Definition of the solar system barycenter (SSB)

The definition of the solar system barycenter at the origin of the time of integration is based on the hypothesis of the conservation of the energy of the dynamical system composed by the planets, the Moon and the asteroids (including Trans-Neptunian Objects) in the solar system. The position *R* and the velocity *V* of the center of mass are then invariant. Using the notation of DE430 (Folkner et al. 2014), it comes

$$ \left\{ \begin{aligned} R &= \frac{\sum_i \mu_i^* \mathbf{r}_i}{\sum_i \mu_i^*} = 0 \\ V &= \frac{dR}{dt} = 0 \end{aligned} \right. \qquad (87) $$

where $\mathbf{r}_i$ and $\mathbf{v}_i$ being respectively its barycentric position and velocity vectors of the planet *i*. $\mu^*$ is given by

$$ \left\{ \begin{aligned} \mu_i^* &= \mu_i \left( 1 + \frac{v_i^2}{2c^2} - \frac{1}{2c^2} \sum_{j \neq i} \frac{\mu_j}{r_{ij}} \right) \\ \dot{\mu}_i^* &= \frac{\mu_i}{2c^2} \left( \sum_{j \neq i} \mu_j \frac{(\mathbf{r}_j - \mathbf{r}_i) \cdot (\mathbf{v}_j + \mathbf{v}_i)}{r_{ij}^3} \right) \end{aligned} \right. \qquad (88) $$

with $\mu_i$ the product of the gravitational constant *G* with the inertial mass of the body *i*. The term $\dot{\mu}_i^*$ is not included in DE430 nor DE440 ephemerides (Folkner et al. 2014; Park et al. 2021) but is accounted for in the INPOP ephemerides (Fienga et al. 2008). The positions and velocities of the SSB are then obtained by integration of the following equations and up to the $c^{-2}$ order,

$$ \left\{ \begin{aligned} \sum_i \mu_i^* \mathbf{r}_i &= 0 \\ \sum_i \mu_i^* \mathbf{v}_i + \dot{\mu}_i^* \mathbf{r}_i &= 0 \end{aligned} \right. \qquad (89) $$

In the INPOP formalism (Fienga et al. 2008), Eqs. (89) and (88) are solved only at the initial step of the planetary integration and a constant vector in positions and velocities is subtracted to all the body positions and velocities for having *R* and *V* to 0. Once the frame is centered on the SSB defined by Eq. (89) at J2000 (initial date of integration for INPOP ephemerides), the equations of motion of the solar system bodies and of the Sun are integrated in this fixed frame. The method described here is
used by the INPOP planetary ephemerides since INPOP08 (Fienga et al. 2008). In JPL DE430 ephemerides (Folkner et al. 2014), the Sun initial coordinates are set up such as *R* and *V* are 0, the equations of motion of the Sun, the Moon, planets, and asteroids being then integrated in the fixed frame. In the former JPL DE versions, such as DE421 (Folkner et al. 2008), the Sun coordinates were set up at each step of integration for maintaining *R* and *V* to 0 all along the integration process. A description of the successive SSB implementations for older JPL versions can be found in Folkner et al. (2014) and discussions concerning possible uncertainties are presented in Fienga et al. (2008) and Folkner et al. (2014). It is interesting to note that, for the tests of the Equivalence Principle (Sect. 4.2), $\mu_i$ of Eqs. (89) and (88) will still correspond to the product of G with the inertial mass of the body. It will then differ from the gravitational mass used for computing the planetary acceleration as in Eq. (90).

## 3.3 The dynamical model

### 3.3.1 Point-mass interactions

In the context of Soffel (2003) and of the metric defined in Eqs. (80) and (84), one can then, in the mass-monopole approximation, write the equations of motion of bodies as well as the conservation laws satisfied by the SSB as given in Damour and Vokrouhlický (1995). In BCRS, the equation of motion for the point-mass interaction is given by

$$ 
\begin{aligned}
\boldsymbol{a}_T = & -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3} \boldsymbol{r}_{AT} \\
& -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3 c^2} \boldsymbol{r}_{AT} \left\{ \begin{aligned}[t]
& v_T^2 + 2v_A^2 - 4\boldsymbol{v}_A \cdot \boldsymbol{v}_T - \frac{3}{2} \left( \frac{\boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A}{r_{AT}} \right)^2 \\
& - \frac{1}{2} \boldsymbol{r}_{AT} \cdot \boldsymbol{a}_A - 4 \sum_{B \neq T} \frac{\mu_B}{r_{TB}} - \sum_{B \neq A} \frac{\mu_B}{r_{AB}} \end{aligned} \right\} \\
& + \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}^3} \left[ 4 \boldsymbol{r}_{AT} \cdot \boldsymbol{v}_T - 3 \boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A \right] (\boldsymbol{v}_T - \boldsymbol{v}_A) \\
& + \frac{7}{2} \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}} \boldsymbol{a}_A,
\end{aligned} 
\quad (90) 
$$

This equation is known as the Einstein–Infeld–Hoffmann–Droste–Lorentz (EIHDL) equation of motion, and it is numerically integrated for obtaining the modern planetary ephemerides in general relativity at the first post-Newtonian $c^{-2}$ level. Additional accelerations induced by the oblateness of the Sun (Sect. 3.3.4) or frame dragging effects (Sect. 3.3.3) are also accounted for in order to describe the motion of the planets at the level of accuracy required by the modern space mission observations. Other multipole moments than the Sun quadrupole moment are negligible at the current level of accuracy.
used by the INPOP planetary ephemerides since INPOP08 (Fienga et al. 2008). In JPL DE430 ephemerides (Folkner et al. 2014), the Sun initial coordinates are set up such as *R* and *V* are 0, the equations of motion of the Sun, the Moon, planets, and asteroids being then integrated in the fixed frame. In the former JPL DE versions, such as DE421 (Folkner et al. 2008), the Sun coordinates were set up at each step of integration for maintaining *R* and *V* to 0 all along the integration process. A description of the successive SSB implementations for older JPL versions can be found in Folkner et al. (2014) and discussions concerning possible uncertainties are presented in Fienga et al. (2008) and Folkner et al. (2014). It is interesting to note that, for the tests of the Equivalence Principle (Sect. 4.2), $\mu_i$ of Eqs. (89) and (88) will still correspond to the product of G with the inertial mass of the body. It will then differ from the gravitational mass used for computing the planetary acceleration as in Eq. (90).

## 3.3 The dynamical model

### 3.3.1 Point-mass interactions

In the context of Soffel (2003) and of the metric defined in Eqs. (80) and (84), one can then, in the mass-monopole approximation, write the equations of motion of bodies as well as the conservation laws satisfied by the SSB as given in Damour and Vokrouhlický (1995). In BCRS, the equation of motion for the point-mass interaction is given by

$$ 
\begin{aligned}
\boldsymbol{a}_T = & -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3} \boldsymbol{r}_{AT} \\
& -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3 c^2} \boldsymbol{r}_{AT} \left\{ \begin{aligned}[t]
& v_T^2 + 2v_A^2 - 4\boldsymbol{v}_A \cdot \boldsymbol{v}_T - \frac{3}{2} \left( \frac{\boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A}{r_{AT}} \right)^2 \\
& - \frac{1}{2} \boldsymbol{r}_{AT} \cdot \boldsymbol{a}_A - 4 \sum_{B \neq T} \frac{\mu_B}{r_{TB}} - \sum_{B \neq A} \frac{\mu_B}{r_{AB}} \end{aligned} \right\} \\
& + \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}^3} \left[ 4\boldsymbol{r}_{AT} \cdot \boldsymbol{v}_T - 3\boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A \right] (\boldsymbol{v}_T - \boldsymbol{v}_A) \\
& + \frac{7}{2} \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}} \boldsymbol{a}_A,
\end{aligned}
\quad (90) 
$$

This equation is known as the Einstein–Infeld–Hoffmann–Droste–Lorentz (EIHDL) equation of motion, and it is numerically integrated for obtaining the modern planetary ephemerides in general relativity at the first post-Newtonian $c^{-2}$ level. Additional accelerations induced by the oblateness of the Sun (Sect. 3.3.4) or frame dragging effects (Sect. 3.3.3) are also accounted for in order to describe the motion of the planets at the level of accuracy required by the modern space mission observations. Other multipole moments than the Sun quadrupole moment are negligible at the current level of accuracy.
### 3.3.2 Asteroid perturbations

Differences in the number of point-mass perturbations exist between planetary ephemerides. Table 2 summarizes them by splitting the case of the Main-Belt asteroids (MBA) on the one hand, and of the Trans-Neptunian objets (TNO) on the other hand. The problem of the modelisation of the MBA perturbations on the inner planet orbits has been addressed since DE405 (Standish 2001). The issue is that there is a considerable number of objects (at least 250,000), for which the masses are unknown, that can potentially interact gravitationally with the inner planets (mainly Mars). With DE405 (Standish 1998a), 5 asteroids (Ceres, Pallas, Vesta, Iris and Bamberga) have been identified as the main perturbers associated with an averaging of the effect induced by the rest of the MBA based on taxonomic classifications. But this modelisation shows its limits with the improvement of the accuracy of the Mars-orbiter tracking data (Standish and Fienga 2002). Williams (1984) proposed a list of 343 individual objects among the MBA gathering the most perturbing asteroids to consider for the construction of accurate planetary ephemerides. That list has been confirmed by Kuchynka and Folkner (2013) and since DE430 this is the Williams (1984) objects that are taken into account in DE and INPOP planetary ephemerides, with their masses individually fitted to observations together with the initial conditions of planet orbits, the mass of the Sun and the Earth-Moon mass ratio (Folkner et al. 2014; Fienga et al. 2019). A detailed description of the complexity of such an adjustment is out of the scope of this review but more comments on the accuracy obtained for these mass determinations can be found in Kuchynka and Folkner (2013) and Fienga et al. (2020a). In order to overcome this problem, but also to reduce the time cost of integrating an important number of objects, it has been proposed to substitute or to complement the MBA individual point-mass contribution by a global ring potential.

The first ring modelisation was proposed by Krasinsky et al. (2002) as a static ring fixed at 2.8 AU from the SSB. An updated model maintaining the total linear and angular momenta of the system has been implemented in Fienga et al. (2009) and Kuchynka et al. (2010) as a way to account for additional MBA contribution associated with a shorter list of individual point-mass perturbers. This approach has not be retained for MBA in DE (Folkner et al. 2014) and INPOP (since INPOP13c (Fienga et al. 2014)) ephemerides but kept in EPM ones associated with individual mass and density determinations for 3 classes of MBA (Pitjeva and Pitjev 2018).

After the introduction of the Juno and Cassini tracking data, it became clear that the perturbations induced by the TNO have to be included in the model. Pitjeva and Pitjev (2018), Di Ruscio et al. (2020b) and Park et al. (2021), individual point-mass perturbations of the most massive and binary Trans-Neptunian objects have been added together with a circular ring for accounting for the average effect of the rest of the TNO. However as the TNO orbits are more eccentric than the MBA ones, Fienga et al. (2021a) had introduced, instead of a circular ring, 500 individual point-mass TNO perturbers with observed orbits<fnref n="19" /> but of equal mass. This latest assumption

<footnote n="19">¹⁹ extracted from the astord database (Moskovitz et al. 2022).</footnote>
### 3.3.2 Asteroid perturbations

Differences in the number of point-mass perturbations exist between planetary ephemerides. Table 2 summarizes them by splitting the case of the Main-Belt asteroids (MBA) on the one hand, and of the Trans-Neptunian objets (TNO) on the other hand. The problem of the modelisation of the MBA perturbations on the inner planet orbits has been addressed since DE405 (Standish 2001). The issue is that there is a considerable number of objects (at least 250,000), for which the masses are unknown, that can potentially interact gravitationally with the inner planets (mainly Mars). With DE405 (Standish 1998a), 5 asteroids (Ceres, Pallas, Vesta, Iris and Bamberga) have been identified as the main perturbers associated with an averaging of the effect induced by the rest of the MBA based on taxonomic classifications. But this modelisation shows its limits with the improvement of the accuracy of the Mars-orbiter tracking data (Standish and Fienga 2002). Williams (1984) proposed a list of 343 individual objects among the MBA gathering the most perturbing asteroids to consider for the construction of accurate planetary ephemerides. That list has been confirmed by Kuchynka and Folkner (2013) and since DE430 this is the Williams (1984) objects that are taken into account in DE and INPOP planetary ephemerides, with their masses individually fitted to observations together with the initial conditions of planet orbits, the mass of the Sun and the Earth-Moon mass ratio (Folkner et al. 2014; Fienga et al. 2019). A detailed description of the complexity of such an adjustment is out of the scope of this review but more comments on the accuracy obtained for these mass determinations can be found in Kuchynka and Folkner (2013) and Fienga et al. (2020a). In order to overcome this problem, but also to reduce the time cost of integrating an important number of objects, it has been proposed to substitute or to complement the MBA individual point-mass contribution by a global ring potential.

The first ring modelisation was proposed by Krasinsky et al. (2002) as a static ring fixed at 2.8 AU from the SSB. An updated model maintaining the total linear and angular momenta of the system has been implemented in Fienga et al. (2009) and Kuchynka et al. (2010) as a way to account for additional MBA contribution associated with a shorter list of individual point-mass perturbers. This approach has not be retained for MBA in DE (Folkner et al. 2014) and INPOP (since INPOP13c (Fienga et al. 2014)) ephemerides but kept in EPM ones associated with individual mass and density determinations for 3 classes of MBA (Pitjeva and Pitjev 2018).

After the introduction of the Juno and Cassini tracking data, it became clear that the perturbations induced by the TNO have to be included in the model. Pitjeva and Pitjev (2018), Di Ruscio et al. (2020b) and Park et al. (2021), individual point-mass perturbations of the most massive and binary Trans-Neptunian objects have been added together with a circular ring for accounting for the average effect of the rest of the TNO. However as the TNO orbits are more eccentric than the MBA ones, Fienga et al. (2021a) had introduced, instead of a circular ring, 500 individual point-mass TNO perturbers with observed orbits<fnref n="19" /> but of equal mass. This latest assumption

<footnote n="19">¹⁹ extracted from the astord database (Moskovitz et al. 2022).</footnote>
simplifies the adjustment as only one parameter (the total mass spread over 500 objects) is used and fitted for characterising the TNO perturbations.

The full 1PN general relativity equation of motion is only taken into account for the first biggest MBA (Ceres, Pallas, Vesta, Iris and Bamberga). For the rest of the MBA and the TNO point-mass interactions, the perturbations on the planet orbits, and on the other asteroid orbits, are estimated at the Newtonian level only. This is done in order to reduce the computational cost, but it is consistent with the overall PN perturbative approach—given that those perturbations are already very small at the Newtonian level, such that their contributions at the post-Newtonian level are negligible.

### 3.3.3 Lense-Thirring acceleration

The Lense-Thirring (LT) acceleration induced by the Sun rotation has been recently added in modern planetary ephemerides (Pitjeva and Pitjev 2018; Park et al. 2021; Fienga et al. 2021b). Its new introduction is justified by the accuracy reached by space missions, especially those close to the Sun and Mercury such as Messenger and BepiColombo (Genova et al. 2018). It also has been shown that it helps for disentangling the PPN parameters for their contribution to the planetary motions (Fienga et al. 2021a).

Iorio et al. (2011) has also shown that this effect contributes to about 10% of the dynamical acceleration induced by the shape of the Sun in General relativity (see Sect. 3.3.4). The acceleration induced by the Lense-Thirring effect generated by a central body (here the Sun) at the first post-Newtonian approximation is given by

$$ \mathbf{a}_{LT} = \frac{2GS}{c^2r^3} \left[ \frac{3k.r}{r^2}(\mathbf{r} \wedge \mathbf{v}) - (\mathbf{k} \wedge \mathbf{v}) \right] \quad (91) $$

where G is the gravitational constant, c the speed of light, S is the Sun angular momentum such as $S = Sk$ where k is the direction of the Sun rotation pole defined according to the IAU right ascension and declination (Archinal et al. 2018), r and v are the position and velocity vectors of the planet relative to the Sun and γ is the PPN parameter for the light deflection, equal to 1 in general relativity. The Lense-Thirring effect hence depends on the Sun angular momentum, S that can be obtained by considering different models for the Sun rotation including both the rotation of the convective region well constrained by helioseismology and the rotation of the Sun core (see Fienga et al. 2021a for discussion).

### 3.3.4 Shape of the sun

Among the solar multipole moments, only the degree 2 order 0 term of the spherical harmonic decomposition of the Sun gravity field, $J_2^\odot$—which is identified to the oblateness of the Sun— leads to a significant impact on the planetary ephemerides at the current level of accuracy. The $J_2^\odot$ induces an acceleration $\mathbf{a}_{J_2^\odot} = (a_{J_2^\odot}^x, a_{J_2^\odot}^y, a_{J_2^\odot}^z)$ to be added to EIHDL Eq. (90) such as (Sharma et al. 2016; Ivashkin 2021):
simplifies the adjustment as only one parameter (the total mass spread over 500 objects) is used and fitted for characterising the TNO perturbations.

The full 1PN general relativity equation of motion is only taken into account for the first biggest MBA (Ceres, Pallas, Vesta, Iris and Bamberga). For the rest of the MBA and the TNO point-mass interactions, the perturbations on the planet orbits, and on the other asteroid orbits, are estimated at the Newtonian level only. This is done in order to reduce the computational cost, but it is consistent with the overall PN perturbative approach—given that those perturbations are already very small at the Newtonian level, such that their contributions at the post-Newtonian level are negligible.

### 3.3.3 Lense-Thirring acceleration

The Lense-Thirring (LT) acceleration induced by the Sun rotation has been recently added in modern planetary ephemerides (Pitjeva and Pitjev 2018; Park et al. 2021; Fienga et al. 2021b). Its new introduction is justified by the accuracy reached by space missions, especially those close to the Sun and Mercury such as Messenger and BepiColombo (Genova et al. 2018). It also has been shown that it helps for disentangling the PPN parameters for their contribution to the planetary motions (Fienga et al. 2021a).

Iorio et al. (2011) has also shown that this effect contributes to about 10% of the dynamical acceleration induced by the shape of the Sun in General relativity (see Sect. 3.3.4). The acceleration induced by the Lense-Thirring effect generated by a central body (here the Sun) at the first post-Newtonian approximation is given by

$$ \mathbf{a}_{LT} = \frac{2GS}{c^2r^3} \left[ \frac{3k.r}{r^2} (\mathbf{r} \wedge \mathbf{v}) - (\mathbf{k} \wedge \mathbf{v}) \right] \quad (91) $$

where G is the gravitational constant, c the speed of light, S is the Sun angular momentum such as $S = Sk$ where k is the direction of the Sun rotation pole defined according to the IAU right ascension and declination (Archinal et al. 2018), r and v are the position and velocity vectors of the planet relative to the Sun and γ is the PPN parameter for the light deflection, equal to 1 in general relativity. The Lense-Thirring effect hence depends on the Sun angular momentum, S that can be obtained by considering different models for the Sun rotation including both the rotation of the convective region well constrained by helioseismology and the rotation of the Sun core (see Fienga et al. 2021a for discussion).

### 3.3.4 Shape of the sun

Among the solar multipole moments, only the degree 2 order 0 term of the spherical harmonic decomposition of the Sun gravity field, $J_2^\odot$—which is identified to the oblateness of the Sun— leads to a significant impact on the planetary ephemerides at the current level of accuracy. The $J_2^\odot$ induces an acceleration $\mathbf{a}_{J_2^\odot} = (a_{J_2^\odot}^x, a_{J_2^\odot}^y, a_{J_2^\odot}^z)$ to be added to EIHDL Eq. (90) such as (Sharma et al. 2016; Ivashkin 2021):
$$ egin{cases} a_{J_2}^x = -\frac{3}{2} \mu_{\odot} J_2^{\odot} R_{\odot}^2 x \frac{(x^2 + y^2 - 4z^2)}{2r^7}, \\ a_{J_2}^y = -\frac{3}{2} \mu_{\odot} J_2^{\odot} R_{\odot}^2 y \frac{(x^2 + y^2 - 4z^2)}{2r^7}, \\ a_{J_2}^z = -\frac{3}{2} \mu_{\odot} J_2^{\odot} R_{\odot}^2 z \frac{(3x^2 + 3y^2 - 2z^2)}{2r^7} \end{cases} \quad (92) $$

where $\mu_{\odot}$ and $R_{\odot}$ are the gravitational parameter and the radius of the Sun, respectively, $r$ is the heliocentric distance of the planet, $(x, y, z)$ being its heliocentric coordinates.

It is common to see in the literature dealing with geodesy and geophysics the contribution of $J_2^{\odot}$ as time variations of the osculating elements of the perturbed orbit such as (Kaula 1966; Bertotti et al. 2003a):

$$ egin{cases} \dot{a} = 0; & \dot{e} = 0; \dot{i} = 0 \\ \dot{\omega} = \frac{3}{4} J_2^{\odot} \frac{n R_{\odot}^2}{(1 - e^2)^2 a^2} (5 \cos^2 i - 1) \\ \dot{\Omega} = -\frac{3}{2} J_2^{\odot} \frac{n R_{\odot}^2}{(1 - e^2)^2 a^2} \cos i \\ \dot{M} = n + \frac{3}{4} J_2^{\odot} \frac{n R_{\odot}^2}{(1 - e^2)^{3/2} a^2} (3 \cos^2 i - 1) \end{cases} \quad (93) $$

with $a, e, n$ and $i$ are the semi-major axis, the eccentricity, the mean motion and the inclination of the planetary orbit impacted by the oblateness of the Sun. $M$ is the planet mean anomaly. We give these equations Eq. (93) as a simple illustration of the type of contribution one may get from the Sun quadrupole. It is then clear that the major contribution of this quadrupole term is on the advance of nodes $\dot{\Omega}$ and perihelia $\dot{\omega}$ of the planetary orbits. This system of equations is also somehow useful to apprehend the potential correlations that may appear between the Sun quadrupole moment and the semi-major axis and eccentricity.

Figure 6 shows $J_2^{\odot}$ values deduced from planetary ephemerides including or not the Lense-Thirring contribution. It is clear, from this figure, that the introduction of the Lense-Thirring acceleration in the model improves significantly the $J_2^{\odot}$ determinations. The values obtained with LT (Fienga et al. 2021a; Park et al. 2017) are indeed closer to the ones issued from helioseismologic surveys (Pijpers 1998; Antia et al. 2008) than those obtained before the LT introduction (Fienga et al. 2019; Iorio et al. 2011; Viswanathan et al. 2017; Fienga et al. 2014; Folkner et al. 2014; Fienga et al. 2011a, 2009). In general relativity, the difference in $J_2^{\odot}$ with or without LT is of about about 7–10%.

As it will be developed in Sect. 4.1, the introduction of gravity theories that are different but close to general relativity induces modifications of the previous Eqs. (90)–(92). In this context, the value of $J_2^{\odot}$ can also be determined with planetary ephemerides in the PPN framework, like given in Table 5. But we will see that the value of $J_2^{\odot}$ can be strongly correlated to the value of non-general relativity
$$
\left\{
\begin{aligned}
a_{J_2}^x &= -\frac{3}{2} \mu_{\odot} J_2^{\odot} R_{\odot}^2 x \frac{(x^2 + y^2 - 4z^2)}{2r^7}, \\
a_{J_2}^y &= -\frac{3}{2} \mu_{\odot} J_2^{\odot} R_{\odot}^2 y \frac{(x^2 + y^2 - 4z^2)}{2r^7}, \\
a_{J_2}^z &= -\frac{3}{2} \mu_{\odot} J_2^{\odot} R_{\odot}^2 z \frac{(3x^2 + 3y^2 - 2z^2)}{2r^7}
\end{aligned}
\right.
\qquad (92)
$$

where $\mu_{\odot}$ and $R_{\odot}$ are the gravitational parameter and the radius of the Sun, respectively, $r$ is the heliocentric distance of the planet, $(x, y, z)$ being its heliocentric coordinates.

It is common to see in the literature dealing with geodesy and geophysics the
contribution of $J_2^{\odot}$ as time variations of the osculating elements of the perturbed orbit
such as (Kaula 1966; Bertotti et al. 2003a):

$$
\left\{
\begin{aligned}
\dot{a} &= 0; & \dot{e} &= 0; & \dot{i} &= 0 \\
\dot{\omega} &= \frac{3}{4} J_2^{\odot} \frac{n}{(1 - e^2)^2} \frac{R_{\odot}^2}{a^2} (5 \cos^2 i - 1) \\
\dot{\Omega} &= -\frac{3}{2} J_2^{\odot} \frac{n}{(1 - e^2)^2} \frac{R_{\odot}^2}{a^2} \cos i \\
\dot{M} &= n + \frac{3}{4} J_2^{\odot} \frac{n}{(1 - e^2)^{3/2}} \frac{R_{\odot}^2}{a^2} (3 \cos^2 i - 1)
\end{aligned}
\right.
\qquad (93)
$$

with *a*, *e*, *n* and *i* are the semi-major axis, the eccentricity, the mean motion and the inclination of the planetary orbit impacted by the oblateness of the Sun. *M* is the planet mean anomaly. We give these equations Eq. (93) as a simple illustration of the type of contribution one may get from the Sun quadrupole. It is then clear that the major contribution of this quadrupole term is on the advance of nodes $\dot{\Omega}$ and perihelia $\dot{\omega}$ of the planetary orbits. This system of equations is also somehow useful to apprehend the potential correlations that may appear between the Sun quadrupole moment and the semi-major axis and eccentricity.

Figure 6 shows $J_2^{\odot}$ values deduced from planetary ephemerides including or not the Lense–Thirring contribution. It is clear, from this figure, that the introduction of the Lense–Thirring acceleration in the model improves significantly the $J_2^{\odot}$ determinations. The values obtained with LT (Fienga et al. 2021a; Park et al. 2017) are indeed closer to the ones issued from helioseismologic surveys (Pijpers 1998; Antia et al. 2008) than those obtained before the LT introduction (Fienga et al. 2019; Iorio et al. 2011; Viswanathan et al. 2017; Fienga et al. 2014; Folkner et al. 2014; Fienga et al. 2011a, 2009). In general relativity, the difference in $J_2^{\odot}$ with or without LT is of about about 7–10%.

As it will be developed in Sect. 4.1, the introduction of gravity theories that are
different but close to general relativity induces modifications of the previous
Eqs. (90)–(92). In this context, the value of $J_2^{\odot}$ can also be determined with planetary
ephemerides in the PPN framework, like given in Table 5. But we will see that the
value of $J_2^{\odot}$ can be strongly correlated to the value of non-general relativity
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Values of the Sun oblateness $J_2$ obtained with different ephemerides and other methods in the general relativity frame. H stands for values deduced from helioseismology, P points values obtained in a partial fit of planetary orbits (i.e. considering only a couple of planets for the fit, like Mercury and the Earth), N indicates values obtained before the introduction of the Lense–Thirring acceleration in the dynamical modeling and Y gives the values obtained after the introduction. The dotted lines give the limits of the less constraining helioseismologic value from Pijpers (1998). In y-axis, are given the references: G18 is Genova et al. (2018), P98 (Pijpers 1998), A08 (Antia et al. 2008), I20a (Fienga et al. 2021a), D440 (Park et al. 2017), I19a (Fienga et al. 2019), I17a (Viswanathan et al. 2017), I13c (Fienga et al. 2014), D414 (Konopliv et al. 2006), I10a (Fienga et al. 2011a) and I08 (Fienga et al. 2009). It is important to note that in Genova et al. (2018), only the orbit of Mercury is considered</figcaption></figure>

parameters, and in particular the PPN parameter $\beta$ (Milani et al. 2002)—as one can get an intuition from Eq. (100). In Sect. 4.1, we will see that the determination of the Sun oblateness when one includes PPN deviations plays a key role in the context of the global adjustment of planetary ephemerides.

## 3.4 Planetary datasets

The modern planetary ephemerides are characterised by the role of planet positions deduced from spacecraft missions either by navigation tracking or by the use of radio science data. In both cases, it is the distance between a station on the ground and a spacecraft orbiting of flying over a planet that is used for the construction of planetary ephemerides. The angular position between the same spacecraft and an extra-galactic source (generally a reference beacon for the construction of the inertial frame, ICRF) is also of great importance for the adjustment of the planetary ephemerides as it allows the tie of the planetary planes to the ICRF (see Sect. 3.2.3).
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Values of the Sun oblateness $J_2$ obtained with different ephemerides and other methods in the general relativity frame. H stands for values deduced from helioseismology, P points values obtained in a partial fit of planetary orbits (i.e. considering only a couple of planets for the fit, like Mercury and the Earth), N indicates values obtained before the introduction of the Lense-Thirring acceleration in the dynamical modeling and Y gives the values obtained after the introduction. The dotted lines give the limits of the less constraining helioseismologic value from Pijpers (1998). In y-axis, are given the references: G18 is Genova et al. (2018), P98 (Pijpers 1998), A08 (Antia et al. 2008), I20a (Fienga et al. 2021a), D440 (Park et al. 2017), I19a (Fienga et al. 2019), I17a (Viswanathan et al. 2017), I13c (Fienga et al. 2014), D414 (Konopliv et al. 2006), I10a (Fienga et al. 2011a) and I08 (Fienga et al. 2009). It is important to note that in Genova et al. (2018), only the orbit of Mercury is considered</figcaption></figure>

parameters, and in particular the PPN parameter $\beta$ (Milani et al. 2002)—as one can get an intuition from Eq. (100). In Sect. 4.1, we will see that the determination of the Sun oblateness when one includes PPN deviations plays a key role in the context of the global adjustment of planetary ephemerides.

## 3.4 Planetary datasets

The modern planetary ephemerides are characterised by the role of planet positions deduced from spacecraft missions either by navigation tracking or by the use of radio science data. In both cases, it is the distance between a station on the ground and a spacecraft orbiting of flying over a planet that is used for the construction of planetary ephemerides. The angular position between the same spacecraft and an extra-galactic source (generally a reference beacon for the construction of the inertial frame, ICRF) is also of great importance for the adjustment of the planetary ephemerides as it allows the tie of the planetary planes to the ICRF (see Sect. 3.2.3).
Optical meridian transit and CCD angular positions are mainly kept for the constraint on the outer planet orbits. A full description of the different data sets used for the construction of the planetary ephemerides can be found in Folkner et al. (2014) and Fienga et al. (2009) and regular updates are presented for every new releases such as for DE440 (Park et al. 2021) or INPOP19a (Di Ruscio et al. 2020b). Table 3 gives the different data sets used for the construction of INPOP21a (Fienga et al. 2021b), counting a total of about 200,000 observations, this number of course varying from one ephemeris to another. Very similar datasets have been used by DE440 (Park et al. 2021).

### 3.4.1 Time delay of the propagation of light

The propagation of light is affected by the curvature of spacetime. Although the trajectory of light is bent in gravitational potentials, this bending is negligible for radio science in the solar system,<fnref n="20" /> and the main effect is known as the Shapiro delay: the propagation time between two points is not simply the Euclidean BCRS distance over the speed of light between an emission point *e* and a reception point *r*, but it also involves a delay in general relativity that depends on the gravitational potentials as follows (Moyer 2000)

$$c(t_r - t_e)_{\text{GRT}} = R + \sum_A 2 \frac{\mu_A}{c^2} \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + \frac{4\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + \frac{4\mu_A}{c^2}}, \quad (94)$$

where the *R* is the Euclidean coordinate distance, $\mathbf{n} = (\mathbf{r}_r - \mathbf{r}_e)/\|\mathbf{r}_r - \mathbf{r}_e\|$, $\mathbf{r}_{iA} = \mathbf{r}_i - \mathbf{r}_A$ and $\mathbf{r}_{eA} = \|\mathbf{r}_{iA}\|$ with $i = e$ or $r$. Note that Eq. (94) slightly differs from the actual Shapiro equation through the $c^{-2}$ terms in the logarithm function. It has lately been realized that this version of the Shapiro delay—rather than the usual equation without the additional $c^{-2}$ terms in the logarithm—allows to account for the so-called second order time propagation enhanced terms, which can become numerically significant for particular conjunction situations despite formally being next-to-leading-order (or $c^{-4}$) terms (Ashby and Bertotti 2010; Linet and Teyssandier 2016). In practice, the difference between the two versions of the Shapiro delay is irrelevant for current planetary ephemerides.

However, it has been realized during the last decade that for close solar conjunctions, additional enhanced $c^{-4}$ terms may no longer be negligible, in particular for BepiColombo MORE (Mercury Orbiter radio science Experiment) (Hees et al. 2014a; Cappuccio et al. 2021) measurements. As a consequence, for close conjunction events, one might have to go beyond Eq. (94) and use a full second order time propagation time formula. This will be challenging in many different ways —see Appendix A. For planetary ephemeris construction, BepiColombo MORE experiment should provide range measurements up to close to 5° to the Sun, that is,

<footnote n="20">²⁰ Because it is a second order effect in the post-Newtonian development of the propagation time.</footnote>
Optical meridian transit and CCD angular positions are mainly kept for the constraint on the outer planet orbits. A full description of the different data sets used for the construction of the planetary ephemerides can be found in Folkner et al. (2014) and Fienga et al. (2009) and regular updates are presented for every new releases such as for DE440 (Park et al. 2021) or INPOP19a (Di Ruscio et al. 2020b). Table 3 gives the different data sets used for the construction of INPOP21a (Fienga et al. 2021b), counting a total of about 200,000 observations, this number of course varying from one ephemeris to another. Very similar datasets have been used by DE440 (Park et al. 2021).

### 3.4.1 Time delay of the propagation of light

The propagation of light is affected by the curvature of spacetime. Although the trajectory of light is bent in gravitational potentials, this bending is negligible for radio science in the solar system,<fnref n="20" /> and the main effect is known as the Shapiro delay: the propagation time between two points is not simply the Euclidean BCRS distance over the speed of light between an emission point *e* and a reception point *r*, but it also involves a delay in general relativity that depends on the gravitational potentials as follows (Moyer 2000)

$$c(t_r - t_e)_{GRT} = R + \sum_A 2 \frac{\mu_A}{c^2} \ln \frac{n \cdot r_{rA} + r_{rA} + \frac{4\mu_A}{c^2}}{n \cdot r_{eA} + r_{eA} + \frac{4\mu_A}{c^2}}, \quad (94)$$

where the *R* is the Euclidean coordinate distance, $n = (r_r - r_e) / \|r_r - r_e\|$, $r_{iA} = r_i - r_A$ and $r_{iA} = \|r_{iA}\|$ with $i = e$ or $r$. Note that Eq. (94) slightly differs from the actual Shapiro equation through the $c^{-2}$ terms in the logarithm function. It has lately been realized that this version of the Shapiro delay—rather than the usual equation without the additional $c^{-2}$ terms in the logarithm—allows to account for the so-called second order time propagation enhanced terms, which can become numerically significant for particular conjunction situations despite formally being next-to-leading-order (or $c^{-4}$) terms (Ashby and Bertotti 2010; Linet and Teyssandier 2016). In practice, the difference between the two versions of the Shapiro delay is irrelevant for current planetary ephemerides.

However, it has been realized during the last decade that for close solar conjunctions, additional enhanced $c^{-4}$ terms may no longer be negligible, in particular for BepiColombo MORE (Mercury Orbiter radio science Experiment) (Hees et al. 2014a; Cappuccio et al. 2021) measurements. As a consequence, for close conjunction events, one might have to go beyond Eq. (94) and use a full second order time propagation time formula. This will be challenging in many different ways —see Appendix A. For planetary ephemeris construction, BepiColombo MORE experiment should provide range measurements up to close to 5° to the Sun, that is,

<footnote n="20">²⁰ Because it is a second order effect in the post-Newtonian development of the propagation time.</footnote>
<table><caption>**Table 3** Detailed example of a planetary ephemeris (here INPOP21a) data sample extracted from Fienga et al. (2021b). Column 1 gives the observed planet and information on the type of observations, Column 2 gives the type of data or the name of spacecraft used for producing the data. Columns 3 and 4 give the time interval and the a priori uncertainties provided by space agencies or the navigation teams, respectively. In Column 2, MRO for Mars Reconnaissance Orbiter and MO for Mars Orbiter. Flybys stands for average normal points obtained during a flyby of a spacecraft with U, C, P, V gives the normal points obtained after the flybys of Ulysses, Cassini, Pioneer and Voyager respectively. radar is for direct radar ranging on planetary surfaces, VLBI stands for Delta-DOR observations and Transit + CCD gives the optical angular observations obtained by transit, photographic plates or CCD. H14 stands for Hees et al. (2014b)</caption><thead><tr><th>Planet/type [unit]</th><th>s/c or method</th><th>Period</th><th>Averaged accuracy</th></tr></thead><tbody><tr><td colspan="4"><strong>Mercury</strong></td></tr><tr><td>Direct range [m]</td><td>Surface</td><td>1971.29: 1997.60</td><td>900</td></tr><tr><td>Radio science range [m]</td><td>Messenger</td><td>2011.23: 2014.26</td><td>5</td></tr><tr><td>Navigation range [m]</td><td>Mariner</td><td>1974.24: 1976.21</td><td>100</td></tr><tr><td colspan="4"><strong>Mercury</strong></td></tr><tr><td>VLBI [mas]</td><td>Magellan, Venus Express</td><td>1990.70: 2013.14</td><td>2.0</td></tr><tr><td>Direct range [m]</td><td>Surface</td><td>1965.96: 1990.07</td><td>1400</td></tr><tr><td>Navigation range [m]</td><td>Venus Express</td><td>2006.32: 2011.45</td><td>7.0</td></tr><tr><td colspan="4"><strong>Mars</strong></td></tr><tr><td>VLBI [mas]</td><td>MGS, MRO</td><td>1989.13: 2013.86</td><td>0.3</td></tr><tr><td>Navigation range [m]</td><td>Mars Express</td><td>2005.17: 2019.37</td><td>2.0</td></tr><tr><td>Radio Science range [m]</td><td>MGS</td><td>1999.31: 2006.70</td><td>2.0</td></tr><tr><td>Radio Science range [m]</td><td>MRO/MO</td><td>2002.14: 2014.00</td><td>1.2</td></tr><tr><td>Navigation range [m]</td><td>Viking</td><td>1976.55: 1982.87</td><td>20.0</td></tr><tr><td colspan="4"><strong>Jupiter</strong></td></tr><tr><td>VLBI [mas]</td><td>Galileo</td><td>1996.54: 1997.94</td><td>11</td></tr><tr><td></td><td>Juno</td><td>2016:2020</td><td>0.5</td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD</td><td>1924.34: 2008.49</td><td>0.3</td></tr><tr><td>Flyby RA/Dec [mas]</td><td>U,C, P, V</td><td>1974.92: 2001.00</td><td>8.0</td></tr><tr><td>Flyby range [m]</td><td>U,C, P, V</td><td>1974.92: 2001.00</td><td>2000</td></tr><tr><td>Radio science range [m]</td><td>Juno</td><td>2016.65: 2020.56</td><td>20</td></tr><tr><td colspan="4"><strong>Saturn</strong></td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD</td><td>1924.22: 2008.34</td><td>0.3</td></tr><tr><td>VLBI RA/Dec [mas]</td><td>Cassini</td><td>2004.69: 2017.9</td><td>0.6</td></tr><tr><td>JPL H14 [m]</td><td>Cassini</td><td>2004.41: 2014.38</td><td>25.0</td></tr><tr><td>Navigation [m]</td><td>Cassini</td><td>2006.01: 2009.83</td><td>6.0</td></tr><tr><td>Radio science: Titan flybys [m]</td><td>Cassini</td><td>2006.01: 2016.61</td><td>15.0</td></tr><tr><td>Radio science: grand finale) [m]</td><td>Cassini</td><td>2017.35: 2017.55</td><td>1.0</td></tr><tr><td colspan="4"><strong>Uranus</strong></td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD+photo</td><td>1924.62: 2013.75</td><td>0.25</td></tr><tr><td>Flyby RA/Dec [mas]</td><td>Voyager</td><td>1986.07</td><td>50</td></tr><tr><td>Flyby range [m]</td><td>Voyager</td><td>1986.07</td><td>50</td></tr><tr><td colspan="4"><strong>Neptune</strong></td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD+photo</td><td>1924.04: 2007.88</td><td>0.3</td></tr><tr><td>Flyby RA/Dec [mas]</td><td>Voyager</td><td>1989.65</td><td>15.0</td></tr><tr><td>Flyby range [m]</td><td>Voyager</td><td>1989.65</td><td>2</td></tr></tbody></table>
<table><caption>**Table 3** Detailed example of a planetary ephemeris (here INPOP21a) data sample extracted from Fienga et al. (2021b). Column 1 gives the observed planet and information on the type of observations, Column 2 gives the type of data or the name of spacecraft used for producing the data. Columns 3 and 4 give the time interval and the a priori uncertainties provided by space agencies or the navigation teams, respectively. In Column 2, MRO for Mars Reconnaissance Orbiter and MO for Mars Orbiter. Flybys stands for average normal points obtained during a flyby of a spacecraft with *U*, *C*, *P*, *V* gives the normal points obtained after the flybys of Ulysses, Cassini, Pioneer and Voyager respectively. radar is for direct radar ranging on planetary surfaces, VLBI stands for Delta-DOR observations and Transit + CCD gives the optical angular observations obtained by transit, photographic plates or CCD. H14 stands for Hees et al. (2014b)</caption><thead><tr><th>Planet/type [unit]</th><th>s/c or method</th><th>Period</th><th>Averaged accuracy</th></tr></thead><tbody><tr><td colspan="4"><strong>Mercury</strong></td></tr><tr><td>Direct range [m]</td><td>Surface</td><td>1971.29: 1997.60</td><td>900</td></tr><tr><td>Radio science range [m]</td><td>Messenger</td><td>2011.23: 2014.26</td><td>5</td></tr><tr><td>Navigation range [m]</td><td>Mariner</td><td>1974.24: 1976.21</td><td>100</td></tr><tr><td colspan="4"><strong>Mercury</strong></td></tr><tr><td>VLBI [mas]</td><td>Magellan, Venus Express</td><td>1990.70: 2013.14</td><td>2.0</td></tr><tr><td>Direct range [m]</td><td>Surface</td><td>1965.96: 1990.07</td><td>1400</td></tr><tr><td>Navigation range [m]</td><td>Venus Express</td><td>2006.32: 2011.45</td><td>7.0</td></tr><tr><td colspan="4"><strong>Mars</strong></td></tr><tr><td>VLBI [mas]</td><td>MGS, MRO</td><td>1989.13: 2013.86</td><td>0.3</td></tr><tr><td>Navigation range [m]</td><td>Mars Express</td><td>2005.17: 2019.37</td><td>2.0</td></tr><tr><td>Radio Science range [m]</td><td>MGS</td><td>1999.31: 2006.70</td><td>2.0</td></tr><tr><td>Radio Science range [m]</td><td>MRO/MO</td><td>2002.14: 2014.00</td><td>1.2</td></tr><tr><td>Navigation range [m]</td><td>Viking</td><td>1976.55: 1982.87</td><td>20.0</td></tr><tr><td colspan="4"><strong>Jupiter</strong></td></tr><tr><td>VLBI [mas]</td><td>Galileo</td><td>1996.54: 1997.94</td><td>11</td></tr><tr><td></td><td>Juno</td><td>2016:2020</td><td>0.5</td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD</td><td>1924.34: 2008.49</td><td>0.3</td></tr><tr><td>Flyby RA/Dec [mas]</td><td>U,C, P, V</td><td>1974.92: 2001.00</td><td>8.0</td></tr><tr><td>Flyby range [m]</td><td>U,C, P, V</td><td>1974.92: 2001.00</td><td>2000</td></tr><tr><td>Radio science range [m]</td><td>Juno</td><td>2016.65: 2020.56</td><td>20</td></tr><tr><td colspan="4"><strong>Saturn</strong></td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD</td><td>1924.22: 2008.34</td><td>0.3</td></tr><tr><td>VLBI RA/Dec [mas]</td><td>Cassini</td><td>2004.69: 2017.9</td><td>0.6</td></tr><tr><td>JPL H14 [m]</td><td>Cassini</td><td>2004.41: 2014.38</td><td>25.0</td></tr><tr><td>Navigation [m]</td><td>Cassini</td><td>2006.01: 2009.83</td><td>6.0</td></tr><tr><td>Radio science: Titan flybys [m]</td><td>Cassini</td><td>2006.01: 2016.61</td><td>15.0</td></tr><tr><td>Radio science: grand finale) [m]</td><td>Cassini</td><td>2017.35: 2017.55</td><td>1.0</td></tr><tr><td colspan="4"><strong>Uranus</strong></td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD+photo</td><td>1924.62: 2013.75</td><td>0.25</td></tr><tr><td>Flyby RA/Dec [mas]</td><td>Voyager</td><td>1986.07</td><td>50</td></tr><tr><td>Flyby range [m]</td><td>Voyager</td><td>1986.07</td><td>50</td></tr><tr><td colspan="4"><strong>Neptune</strong></td></tr><tr><td>Optical RA/Dec [arcsec]</td><td>Transit+CCD+photo</td><td>1924.04: 2007.88</td><td>0.3</td></tr><tr><td>Flyby RA/Dec [mas]</td><td>Voyager</td><td>1989.65</td><td>15.0</td></tr><tr><td>Flyby range [m]</td><td>Voyager</td><td>1989.65</td><td>2</td></tr></tbody></table>
close to conjunction (Fienga et al. 2021a). For such close-to-the-Sun observations, a specific analysis will have to be done in order to account for the neglected terms of Eq. (94), but in the framework defined by the IAU (harmonic gauge).

In any case, when one deals with an alternative theory of gravity, one has to take care of the consistency between the different levels of construction of the ephemeris: modifying the Shapiro delay or the equation of motion for massive bodies is not without consequences on the other.

### 3.4.2 Radio science and navigation data

This type of data represent more than 65% of the full data sample used for the construction of the planetary ephemerides. It gathers measurements of distances in terms of time delays obtained during orbital or flyby phases of space missions, either during the regular navigation process, or during specific scientific sessions dedicated to radio science experiments. In this latter case, a detailed assessment of the capabilities of the radio transpondeur used for the re-emission of the captured signal is needed for the analysis (see e.g. Cappuccio et al. 2020). The measurement is the time differences between the time of emission of the signal by a ground-based antenna, the time of capture by the transponder on board of a spacecraft, the time of re-emission by this transponder after amplification and the time of reception of the transmitted signal by the same or a different ground-based antenna. Using a known orbit of the spacecraft around or near the planet, it is then possible to deduce the planetary range measurement, by accounting for the distance between the spacecraft and the centre of mass of the planet estimated during the s/c orbit determination. The accuracy of the range constraints deduced from the radio science and navigation data is then strongly dependent on the uncertainties of the spacecraft orbit determination.

The distribution of the data per missions is given in Fig. 7. In numbers, the Mars orbiter observations dominate, followed by the Venus tracking data. In terms of accuracy, the Mars data also lead the fit together with the less numerous but important tracking data of the Juno mission orbiting Jupiter as well as the Cassini radio science and tracking data for Saturn. The predominance in numbers and accuracy of the Mars observations on the planetary ephemerides data sample explain why some tests of general relativity have been said to be driven by the Mars orbit (Will 2018a), whereas Bernus et al. (2019) stressed also the importance of the outer planet radio science experiments.

Furthermore, in terms of gravitational framework, these observations are also given in the framework used for studying the motion of the spacecraft around the planet. This framework is usually general relativity and discussions have been held regarding the consistency of testing alternative theories of gravity with observational constraints obtained in general relativity (see Verma et al. 2014). The most simple way to solve this problem is to operate the spacecraft orbit determination (OD) in the
close to conjunction (Fienga et al. 2021a). For such close-to-the-Sun observations, a specific analysis will have to be done in order to account for the neglected terms of Eq. (94), but in the framework defined by the IAU (harmonic gauge).

In any case, when one deals with an alternative theory of gravity, one has to take care of the consistency between the different levels of construction of the ephemeris: modifying the Shapiro delay or the equation of motion for massive bodies is not without consequences on the other.

### 3.4.2 Radio science and navigation data

This type of data represent more than 65% of the full data sample used for the construction of the planetary ephemerides. It gathers measurements of distances in terms of time delays obtained during orbital or flyby phases of space missions, either during the regular navigation process, or during specific scientific sessions dedicated to radio science experiments. In this latter case, a detailed assessment of the capabilities of the radio transpondeur used for the re-emission of the captured signal is needed for the analysis (see e.g. Cappuccio et al. 2020). The measurement is the time differences between the time of emission of the signal by a ground-based antenna, the time of capture by the transponder on board of a spacecraft, the time of re-emission by this transponder after amplification and the time of reception of the transmitted signal by the same or a different ground-based antenna. Using a known orbit of the spacecraft around or near the planet, it is then possible to deduce the planetary range measurement, by accounting for the distance between the spacecraft and the centre of mass of the planet estimated during the s/c orbit determination. The accuracy of the range constraints deduced from the radio science and navigation data is then strongly dependent on the uncertainties of the spacecraft orbit determination.

The distribution of the data per missions is given in Fig. 7. In numbers, the Mars orbiter observations dominate, followed by the Venus tracking data. In terms of accuracy, the Mars data also lead the fit together with the less numerous but important tracking data of the Juno mission orbiting Jupiter as well as the Cassini radio science and tracking data for Saturn. The predominance in numbers and accuracy of the Mars observations on the planetary ephemerides data sample explain why some tests of general relativity have been said to be driven by the Mars orbit (Will 2018a), whereas Bernus et al. (2019) stressed also the importance of the outer planet radio science experiments.

Furthermore, in terms of gravitational framework, these observations are also given in the framework used for studying the motion of the spacecraft around the planet. This framework is usually general relativity and discussions have been held regarding the consistency of testing alternative theories of gravity with observational constraints obtained in general relativity (see Verma et al. 2014). The most simple way to solve this problem is to operate the spacecraft orbit determination (OD) in the
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Example of distribution in percentages of the radio science and tracking data and the optical observations. On the left-hand side, data samples are identified by planets (Ma stand for Mars, Me for Mercury, Ve for Venus, Ju for Jupiter, Sa Saturn and UN for Uranus and Neptune) and by missions. On the right-hand side, the data sample are identified according to the observed planet. A specific tag is made on the set of astrometric positions obtained by Camargo et al. (2022). The presented datasets correspond to the one of INPOP21a (see Table 3). Roughly similar distribution is present for DE and EPM ephemerides</figcaption></figure>

same framework as the one of the alternative theories to be tested. But as the s/c navigation is regularly interrupted, even during orbital phases, by engine boosts or attitude corrections, introducing significant changes in the orbit, OD is usually done over a short duration of time (from few hours to few days), in between these manoeuvres. Most of the time, up to now, alternative theories of gravity start to significantly differ from general relativity over longer time intervals. But it is clear that with the increase of the accuracy in the s/c orbit tracking, the modified general relativity that one wants to test should also be implemented during the OD process. See for example De Marchi and Cascioli (2020).

In parallel to the radar tracking, navigation teams also acquire positions of the spacecraft relative to extra-galactic sources by using facilities also dedicated to VLBI observations. This type of observations, called Delta-DOR, represent less than 1% of the full data sample but they are crucial for linking the ephemeris reference frame to the ICRF (see Sect. 3.2.3). The first were obtained in 1981 (Donivan and Newhall 1981) with the Viking Mars lander. Technical description related to the production of this type of observations can be found in Folkner et al. (1994) and Moyer (2000). The most recent planetary ephemerides account for a total of about 300 VLBI positions, mainly extracted from the tracking data of Mars, Venus, Jupiter and Saturn orbiters with accuracies about less than 1 milliarcsecond of degree (mas) (Jones et al. 2011, 2019; Park et al. 2021).
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Example of distribution in percentages of the radio science and tracking data and the optical observations. On the left-hand side, data samples are identified by planets (Ma stand for Mars, Me for Mercury, Ve for Venus, Ju for Jupiter, Sa Saturn and UN for Uranus and Neptune) and by missions. On the right-hand side, the data sample are identified according to the observed planet. A specific tag is made on the set of astrometric positions obtained by Camargo et al. (2022). The presented datasets correspond to the one of INPOP21a (see Table 3). Roughly similar distribution is present for DE and EPM ephemerides</figcaption></figure>

same framework as the one of the alternative theories to be tested. But as the s/c navigation is regularly interrupted, even during orbital phases, by engine boosts or attitude corrections, introducing significant changes in the orbit, OD is usually done over a short duration of time (from few hours to few days), in between these manoeuvres. Most of the time, up to now, alternative theories of gravity start to significantly differ from general relativity over longer time intervals. But it is clear that with the increase of the accuracy in the s/c orbit tracking, the modified general relativity that one wants to test should also be implemented during the OD process. See for example De Marchi and Cascioli (2020).

In parallel to the radar tracking, navigation teams also acquire positions of the spacecraft relative to extra-galactic sources by using facilities also dedicated to VLBI observations. This type of observations, called Delta-DOR, represent less than 1% of the full data sample but they are crucial for linking the ephemeris reference frame to the ICRF (see Sect. 3.2.3). The first were obtained in 1981 (Donivan and Newhall 1981) with the Viking Mars lander. Technical description related to the production of this type of observations can be found in Folkner et al. (1994) and Moyer (2000). The most recent planetary ephemerides account for a total of about 300 VLBI positions, mainly extracted from the tracking data of Mars, Venus, Jupiter and Saturn orbiters with accuracies about less than 1 milliarcsecond of degree (mas) (Jones et al. 2011, 2019; Park et al. 2021).
### 3.4.3 Angular ground-based observations

Finally, the observational data base is completed by the optical observations of outer planets. This set constitutes about 35% of the full data sample, and, despite their uncertainties, are crucial for Uranus and Neptune. For these two planets, only 1 flyby per planet has been acquired during the Voyager missions, leading to a weak constraint for several decade-long orbits. Direct transit observations of planets or satellites are part of the optical data samples as well as astrometric positions deduced from CCD observations and photographic plates. Both planet and satellite positions are included. In the case of satellite astrometry, the position of the planetary barycenter is deduced using satellite ephemerides (see Camargo et al. 2022). The long-term (several decades) behaviour of outer planet orbits is then driven by the planet barycenter equation of motion as given in Eq. (96), but also by the satellite dynamics. Additionally, depending on the reference frames used for the data publications, different algorithms are applied for transforming the astrometric positions (for example given in FK3, FK4 or FK5) to ICRF. Such transformations guarantee, at the level of accuracy of the optical data (about 50 to 100 mas), the link between the INPOP outer planet frame and ICRF. VLBI observations enforce this link at the VLBI accuracy which means a factor of at least 100 compared to the optical tie. On Fig. 7 are given an example (for INPOP21a) of the percentages of optical observations per planet. A specific tag is made on the data sample provided by Camargo et al. (2022) who did a new reduction of astrometric long-term observations of the Uranian main satellites obtained at the Pico dos Dias from 1982 to 2011, using the Gaia EDR3 as reference.

## 3.5 Fitting procedures

The main construction of planetary ephemerides uses a classic least squares approach. The total number of fitted parameters can be found in Table 2. It changes from one ephemeris to another, depending the perturbations included in the model and the observational parametrizations required for the data analysis. Thus, to the parameters presented in Table 2, spacecraft bias accounting for additional transponder delays or additional ground-station delays (Kuchynka et al. 2012) have also to be added as well as parameters related to solar plasma corrections (Verma et al. 2013; Pitjeva et al. 2021).

Complementary methods have been proposed to make accuracy assessments, to overcome the problems of high correlated parameters, in particular in the frame of testing alternative theories of gravity (see Sect. 4), or of multiple correlations between parameters, as in the case of the determination of asteroid masses (Fienga et al. 2020a). In the context of computing threshold values for the violation of general relativity, Fienga et al. (2015) had tested genetic algorithm approaches for identifying sets of parameters (PPN $\beta, \gamma$, the Sun oblateness $J_2^*$ and secular variations of the gravitational mass of the Sun, $\dot{\mu}/\mu$) with which planetary ephemerides can be computed and fitted to the observations with a comparable accuracy to the ephemeris (in this case INPOP15a) built in general relativity. Different approaches
### 3.4.3 Angular ground-based observations

Finally, the observational data base is completed by the optical observations of outer planets. This set constitutes about 35% of the full data sample, and, despite their uncertainties, are crucial for Uranus and Neptune. For these two planets, only 1 flyby per planet has been acquired during the Voyager missions, leading to a weak constraint for several decade-long orbits. Direct transit observations of planets or satellites are part of the optical data samples as well as astrometric positions deduced from CCD observations and photographic plates. Both planet and satellite positions are included. In the case of satellite astrometry, the position of the planetary barycenter is deduced using satellite ephemerides (see Camargo et al. 2022). The long-term (several decades) behaviour of outer planet orbits is then driven by the planet barycenter equation of motion as given in Eq. (96), but also by the satellite dynamics. Additionally, depending on the reference frames used for the data publications, different algorithms are applied for transforming the astrometric positions (for example given in FK3, FK4 or FK5) to ICRF. Such transformations guarantee, at the level of accuracy of the optical data (about 50 to 100 mas), the link between the INPOP outer planet frame and ICRF. VLBI observations enforce this link at the VLBI accuracy which means a factor of at least 100 compared to the optical tie. On Fig. 7 are given an example (for INPOP21a) of the percentages of optical observations per planet. A specific tag is made on the data sample provided by Camargo et al. (2022) who did a new reduction of astrometric long-term observations of the Uranian main satellites obtained at the Pico dos Dias from 1982 to 2011, using the Gaia EDR3 as reference.

## 3.5 Fitting procedures

The main construction of planetary ephemerides uses a classic least squares approach. The total number of fitted parameters can be found in Table 2. It changes from one ephemeris to another, depending the perturbations included in the model and the observational parametrizations required for the data analysis. Thus, to the parameters presented in Table 2, spacecraft bias accounting for additional transponder delays or additional ground-station delays (Kuchynka et al. 2012) have also to be added as well as parameters related to solar plasma corrections (Verma et al. 2013; Pitjeva et al. 2021).

Complementary methods have been proposed to make accuracy assessments, to overcome the problems of high correlated parameters, in particular in the frame of testing alternative theories of gravity (see Sect. 4), or of multiple correlations between parameters, as in the case of the determination of asteroid masses (Fienga et al. 2020a). In the context of computing threshold values for the violation of general relativity, Fienga et al. (2015) had tested genetic algorithm approaches for identifying sets of parameters (PPN $β, γ$, the Sun oblateness $J_2^*$ and secular variations of the gravitational mass of the Sun, $\dot{\mu}/\mu$) with which planetary ephemerides can be computed and fitted to the observations with a comparable accuracy to the ephemeris (in this case INPOP15a) built in general relativity. Different approaches
based on $\chi^2$ analysis with fixed PPN or alternative theory parameters have been investigated as well, using random walk exploration methods associated with various cost functions (Bernus et al. 2019, 2020, 2022; Fienga et al. 2021a). These methods have in common to be very conservative and to give larger constraints than those obtained with the least square procedures. See Sect. 4.1 for the full discussion.

*A final but important comment is about the masses used for the construction of the ephemerides*

Most of the planetary masses are not obtained during the construction of the planetary ephemerides but during the orbit determination of one or several spacecraft flying in planetary systems (see Table 4). They are obtained by combination of data from different missions and techniques, usually in the general relativity framework. For example, in the case of Neptune and Uranus, it is even a combination of Voyager radio tracking data, ground-based optical observations and images acquired by the s/c that leads to the determination of the mass of the planetary system. The s/c orbit and gravity field (mass) determinations are complex as they depend on various sources of interactions from solar radiation pressure, atmospheric dragging, maneuver residual accelerations to non-stochastic accelerations. For that reason, these s/c orbit and gravity field determinations are usually done independently from the planetary ephemerides.<fnref n="21" /> These masses are therefore obtained in the framework of the s/c orbit and of the gravity field determination algorithm, which assumes generally general relativity. Only, the gravitational mass of the Sun, the ratio between the mass of the Earth and the mass of the Moon (invariant to any change of the gravitational constant) as well as a variable number of asteroid gravitational masses (see Sect. 3.3.2) are, in fact, estimated during the planetary ephemeris adjustment. In the case of a test of the violation of the Equivalence principle, the fact that some of the masses used in the ephemeris are obtained in the general relativity framework can introduce an inconsistency in the method. In principle, it should be necessary to consider a more global approach, including the s/c orbit determination in the same framework than the one of the planetary ephemeris. So far, such an attempt was only operated for Parametrized post-Newtonian tests (see Sect. 4.1).

**4 Tests of alternative theoretical frameworks with planetary ephemeris**

Thanks to the accuracy of the most recent observations but also to the more than one century long data span, planetary ephemerides are excellent tools for testing general relativity and its alternatives. In principle, if the model of gravity is not accurate enough, or does not well represent the physics as seen by the observations, it should then lead to a deterioration of the differences (residuals) between observed and computed quantities—deduced from the ephemeris. When the deterioration of residuals built with an alternative theory becomes statistically significant with respect to the residuals obtained in the general relativity framework, one can say that it disfavors the alternative theoretical framework—or favors general relativity over its

<footnote n="21">²¹ Even if iterations are sometimes necessary.</footnote>
based on $\chi^2$ analysis with fixed PPN or alternative theory parameters have been investigated as well, using random walk exploration methods associated with various cost functions (Bernus et al. 2019, 2020, 2022; Fienga et al. 2021a). These methods have in common to be very conservative and to give larger constraints than those obtained with the least square procedures. See Sect. 4.1 for the full discussion.

*A final but important comment is about the masses used for the construction of the ephemerides*

Most of the planetary masses are not obtained during the construction of the planetary ephemerides but during the orbit determination of one or several spacecraft flying in planetary systems (see Table 4). They are obtained by combination of data from different missions and techniques, usually in the general relativity framework. For example, in the case of Neptune and Uranus, it is even a combination of Voyager radio tracking data, ground-based optical observations and images acquired by the s/c that leads to the determination of the mass of the planetary system. The s/c orbit and gravity field (mass) determinations are complex as they depend on various sources of interactions from solar radiation pressure, atmospheric dragging, maneuver residual accelerations to non-stochastic accelerations. For that reason, these s/c orbit and gravity field determinations are usually done independently from the planetary ephemerides.<fnref n="21" /> These masses are therefore obtained in the framework of the s/c orbit and of the gravity field determination algorithm, which assumes generally general relativity. Only, the gravitational mass of the Sun, the ratio between the mass of the Earth and the mass of the Moon (invariant to any change of the gravitational constant) as well as a variable number of asteroid gravitational masses (see Sect. 3.3.2) are, in fact, estimated during the planetary ephemeris adjustment. In the case of a test of the violation of the Equivalence principle, the fact that some of the masses used in the ephemeris are obtained in the general relativity framework can introduce an inconsistency in the method. In principle, it should be necessary to consider a more global approach, including the s/c orbit determination in the same framework than the one of the planetary ephemeris. So far, such an attempt was only operated for Parametrized post-Newtonian tests (see Sect. 4.1).

## 4 Tests of alternative theoretical frameworks with planetary ephemeris

Thanks to the accuracy of the most recent observations but also to the more than one century long data span, planetary ephemerides are excellent tools for testing general relativity and its alternatives. In principle, if the model of gravity is not accurate enough, or does not well represent the physics as seen by the observations, it should then lead to a deterioration of the differences (residuals) between observed and computed quantities—deduced from the ephemeris. When the deterioration of residuals built with an alternative theory becomes statistically significant with respect to the residuals obtained in the general relativity framework, one can say that it disfavors the alternative theoretical framework—or favors general relativity over its

<footnote n="21">²¹ Even if iterations are sometimes necessary.</footnote>
alternative. On the opposite, if the residuals obtained with an alternative theory are significantly improved—that is, if there are significantly smaller—in comparison to the one obtained in the general relativity framework, it favors the alternative theory. However, so far, this has never occurred.

In the following sections, we will explore various theories or phenomenologies that have been scrutinized using planetary ephemerides. The selection of certain theories or phenomenologies to test, rather than others, can be driven by several considerations. One of them could be the historical circumstances. For example, the Parametrized Post-Newtonian (PPN) framework was developed during the early stages of the general relativity testing era and has since been widely employed in tests involving planetary ephemerides.

Another factor could be the complexity of the task, both in terms of numerical computations and statistical methodology. The larger the number of parameters to constrain a theory has, the larger the dimension of the theoretical space there is to examine. Additionally, there is an increasing risk of introducing multiple correlations between the parameters of the theory, as well as between these parameters and those of the solar system (i.e., masses). Given that it can take up to eight hours of computation to adjust the solar system parameters of the ephemeris at a single point in the parameter space (Mariani et al. 2023), achieving a densely filled parameter space becomes highly computationally demanding if the parameter space dimension is greater than one. Moreover, the statistical methodology also grows more complex with an increase in the number of parameters to constrain.

There is also the risk of missing a genuine detection—that is, obtaining better residuals in an alternative theory than in general relativity for a specific parameter value—if the parameter space of alternative theories explored is too sparsely populated. Therefore, it is of paramount importance to have a densely filled parameter space, which consequently makes computational time requirements escalate with the number of parameters to constrain. In the case of a classic least square approach, the risk of high correlations between parameters associated with the observational uncertainties also limit the possible theories to be tested. As a result, tests are often conducted on theories or phenomenologies that hinge on few parameters only.

## 4.1 Parametrized post-Newtonian framework

For many metric theories—that is, theories for which the additional gravitational fields<fnref n="22" /> do not couple to matter directly, such that the theories satisfy the weak equivalence principle by design—all the differences with respect to general relativity can be taken into account through parameters that appear in front of the gravitational potentials in the metric (Will 2014, 2018b). This is the case for instance for scalar-tensor theories, vector-tensor theories, or tensor-vector-scalar theories, and notably TeVeS (Will 2014). Indeed, because the weak equivalence principle is enforced by hand by imposing that the additional gravitational fields do not couple to matter

<footnote n="22">²² Sometimes called “degrees of freedom” in the theoretical literature.</footnote>
alternative. On the opposite, if the residuals obtained with an alternative theory are significantly improved—that is, if there are significantly smaller—in comparison to the one obtained in the general relativity framework, it favors the alternative theory. However, so far, this has never occurred.

In the following sections, we will explore various theories or phenomenologies that have been scrutinized using planetary ephemerides. The selection of certain theories or phenomenologies to test, rather than others, can be driven by several considerations. One of them could be the historical circumstances. For example, the Parametrized Post-Newtonian (PPN) framework was developed during the early stages of the general relativity testing era and has since been widely employed in tests involving planetary ephemerides.

Another factor could be the complexity of the task, both in terms of numerical computations and statistical methodology. The larger the number of parameters to constrain a theory has, the larger the dimension of the theoretical space there is to examine. Additionally, there is an increasing risk of introducing multiple correlations between the parameters of the theory, as well as between these parameters and those of the solar system (i.e., masses). Given that it can take up to eight hours of computation to adjust the solar system parameters of the ephemeris at a single point in the parameter space (Mariani et al. 2023), achieving a densely filled parameter space becomes highly computationally demanding if the parameter space dimension is greater than one. Moreover, the statistical methodology also grows more complex with an increase in the number of parameters to constrain.

There is also the risk of missing a genuine detection—that is, obtaining better residuals in an alternative theory than in general relativity for a specific parameter value—if the parameter space of alternative theories explored is too sparsely populated. Therefore, it is of paramount importance to have a densely filled parameter space, which consequently makes computational time requirements escalate with the number of parameters to constrain. In the case of a classic least square approach, the risk of high correlations between parameters associated with the observational uncertainties also limit the possible theories to be tested. As a result, tests are often conducted on theories or phenomenologies that hinge on few parameters only.

## 4.1 Parametrized post-Newtonian framework

For many metric theories—that is, theories for which the additional gravitational fields<fnref n="22" /> do not couple to matter directly, such that the theories satisfy the weak equivalence principle by design—all the differences with respect to general relativity can be taken into account through parameters that appear in front of the gravitational potentials in the metric (Will 2014, 2018b). This is the case for instance for scalar-tensor theories, vector-tensor theories, or tensor-vector-scalar theories, and notably TeVeS (Will 2014). Indeed, because the weak equivalence principle is enforced by hand by imposing that the additional gravitational fields do not couple to matter

<footnote n="22">²² Sometimes called “degrees of freedom” in the theoretical literature.</footnote>
<table><caption>**Table 4** Planetary masses used in DE440 (Park et al. 2021) and INPOP19a (Fienga et al. 2021b): List of sources providing masses used in the ephemerides, with Column 1 giving the considered planet, Column 2 the reference and Column 3 the type of observations used for the mass computation</caption><thead><tr><th>Planet</th><th>References</th><th>Mission</th></tr></thead><tbody><tr><td>Mercury</td><td>Konopliv et al. (2020)</td><td>Messenger</td></tr><tr><td>Venus</td><td>Konopliv et al. (1999)</td><td>MAGELLAN</td></tr><tr><td>Mars</td><td>Konopliv et al. (2016)</td><td>MRO, MGS, Odyssey, Pathfinder, MER, Viking</td></tr><tr><td>Jupiter barycenter</td><td>Durante et al. (2020)<br>Jacobson (2021b)</td><td>Juno + 5 s/c flybys + ground-based and s/c optical data</td></tr><tr><td>Saturn barycenter</td><td>Jacobson (2021a)<br>Di Ruscio et al. (2020a)</td><td>Cassini + ground-based + s/c optical data</td></tr><tr><td>Uranus barycenter</td><td>Jacobson (2014)</td><td>Voyager 2 + ground-based optical data</td></tr><tr><td>Neptune barycenter</td><td>Jacobson (2009)</td><td>Voyager + ground-based and s/c optical data</td></tr></tbody></table>

directly, particles still follow geodesics of the metric and proper time is still defined from the metric only.<fnref n="23" />

Many theories then turn out to generate the same type of metrics, where the differences with respect to general relativity can be parametrized by a set of a few parameters in the metric only—which values depend upon the theory considered. This comes in handy, as it allows to derive all the equations necessary in planetary ephemerides—that are, the time-scale definitions, the motion of celestial bodies and the Shapiro delay—from a single original parametrized metric. This whole framework is therefore called the parametrized post-Newtonian (PPN) formalism.

Nevertheless, as it will be discussed in Sect. 4.1.2, the coordinate system being used for planetary ephemerides matters, and one has to make sure that the coordinate system being used for alternative theories is compatible at the required level of accuracy with the IAU standards used to convert observables (e.g. the roundtrip of an electromagnetic signal measured in terms of the proper time of the local clock) into spatial and temporal coordinate positions (e.g. the distance<fnref n="24" /> of a spacecraft with respect to the ground-based station being used).

Among the 10 PPN parameters discussed in Will (2014) and Will (2018b), we shall first only consider the Eddington–Robertson–Schiff parameters $\gamma$ and $\beta$, because most of the planetary ephemerides to this date focused on those parameters. A complete description of the PPN parameters will be discussed in Sect. 4.1.4. Let us write the metric as follows (Klioner and Soffel 2000)

<footnote n="23">²³ Technically, a theory satisfies the weak equivalence principle if it exists a conformal (or disformal) representation for which the metric is the only gravitational field that appears in the matter Lagrangian density according to the *comma-goes-to-semicolon rule* (Misner et al. 1973). For instance, the additional scalar-field appears in the matter Lagrangian density in the Brans–Dicke theory when written in the Einstein representation (or frame), but not in the (original) Jordan representation. Hence, because of the latter, the Brans–Dicke theory satisfies the weak equivalence principle.</footnote>

<footnote n="24">²⁴ Which is dependent upon the coordinate system being used.</footnote>
<table><caption>**Table 4** Planetary masses used in DE440 (Park et al. 2021) and INPOP19a (Fienga et al. 2021b): List of sources providing masses used in the ephemerides, with Column 1 giving the considered planet, Column 2 the reference and Column 3 the type of observations used for the mass computation</caption><thead><tr><th>Planet</th><th>References</th><th>Mission</th></tr></thead><tbody><tr><td>Mercury</td><td>Konopliv et al. (2020)</td><td>Messenger</td></tr><tr><td>Venus</td><td>Konopliv et al. (1999)</td><td>MAGELLAN</td></tr><tr><td>Mars</td><td>Konopliv et al. (2016)</td><td>MRO, MGS, Odyssey, Pathfinder, MER, Viking</td></tr><tr><td>Jupiter barycenter</td><td>Durante et al. (2020)<br>Jacobson (2021b)</td><td>Juno + 5 s/c flybys + ground-based and s/c optical data</td></tr><tr><td>Saturn barycenter</td><td>Jacobson (2021a)<br>Di Ruscio et al. (2020a)</td><td>Cassini + ground-based + s/c optical data</td></tr><tr><td>Uranus barycenter</td><td>Jacobson (2014)</td><td>Voyager 2 + ground-based optical data</td></tr><tr><td>Neptune barycenter</td><td>Jacobson (2009)</td><td>Voyager + ground-based and s/c optical data</td></tr></tbody></table>

directly, particles still follow geodesics of the metric and proper time is still defined from the metric only.<fnref n="23" />

Many theories then turn out to generate the same type of metrics, where the differences with respect to general relativity can be parametrized by a set of a few parameters in the metric only—which values depend upon the theory considered. This comes in handy, as it allows to derive all the equations necessary in planetary ephemerides—that are, the time-scale definitions, the motion of celestial bodies and the Shapiro delay—from a single original parametrized metric. This whole framework is therefore called the parametrized post-Newtonian (PPN) formalism.

Nevertheless, as it will be discussed in Sect. 4.1.2, the coordinate system being used for planetary ephemerides matters, and one has to make sure that the coordinate system being used for alternative theories is compatible at the required level of accuracy with the IAU standards used to convert observables (e.g. the roundtrip of an electromagnetic signal measured in terms of the proper time of the local clock) into spatial and temporal coordinate positions (e.g. the distance<fnref n="24" /> of a spacecraft with respect to the ground-based station being used).

Among the 10 PPN parameters discussed in Will (2014) and Will (2018b), we shall first only consider the Eddington–Robertson–Schiff parameters $\gamma$ and $\beta$, because most of the planetary ephemerides to this date focused on those parameters. A complete description of the PPN parameters will be discussed in Sect. 4.1.4. Let us write the metric as follows (Klioner and Soffel 2000)

<footnote n="23">²³ Technically, a theory satisfies the weak equivalence principle if it exists a conformal (or disformal) representation for which the metric is the only gravitational field that appears in the matter Lagrangian density according to the *comma-goes-to-semicolon rule* (Misner et al. 1973). For instance, the additional scalar-field appears in the matter Lagrangian density in the Brans–Dicke theory when written in the Einstein representation (or frame), but not in the (original) Jordan representation. Hence, because of the latter, the Brans–Dicke theory satisfies the weak equivalence principle.</footnote>

<footnote n="24">²⁴ Which is dependent upon the coordinate system being used.</footnote>
$$ 
\begin{aligned}
g_{00} &= -1 + \frac{2w}{c^2} - \frac{2\beta w^2}{c^4} + \mathcal{O}(c^{-5}) \\
g_{0i} &= -\frac{(4\gamma + 3)w_i}{2c^3} + \mathcal{O}(c^{-4}) \\
g_{ij} &= \delta_{ij} \left(1 + \frac{2\gamma w}{c^2}\right) + \mathcal{O}(c^{-3}),
\end{aligned} 
\quad (95) $$

In general relativity, the PPN parameters $\beta$ and $\gamma$ equal to one, but not for alternative theories in general.

The form of the metric given in Eq. (95) is not sufficient in order to completely fix the coordinate system, and one has to impose an additional constraint (or gauge) in order to specify the metric field equations—see, e.g., Klioner and Soffel (2000). Since the IAU recommends the coordinates (or gauge) to be harmonic for general relativity (Soffel 2003), it seems relevant to consider harmonic coordinates for alternative theories as well, such that (at least) the coordinate systems of the alternative theory reduce to the IAU coordinate systems in the limit for which both theories cannot be distinguished at the level of accuracy of the observations.

### Ambiguity of the definition of harmonic coordinates

The very definition of what one may call harmonic coordinates can be ambiguous for alternative theories of gravity. For instance, for scalar-tensor theories, one can impose the usual post-Newtonian harmonic conditions—$g^{\alpha\beta}\Gamma^0_{\alpha\beta} = 0$—either in the Jordan frame, or in the conformal Einstein frame (Damour and Esposito-Farese 1992; Kopeikin and Vlasov 2004; Minazzoli and Chauvineau 2011; Will 2018b; Kopeikin 2019). But both choices reduce to the general relativity harmonic coordinate in the limit for which the difference between general relativity and the scalar-tensor theory becomes negligible (e.g. for the Brans–Dicke parameter that is such that $\omega \to \infty$). Imposing the harmonic conditions in the Einstein frame has been shown to lead to more convenient properties (Damour and Esposito-Farese 1992; Kopeikin and Vlasov 2004; Minazzoli and Chauvineau 2011; Kopeikin 2019), notably with respect to the full theory of reference frames (Kopeikin and Vlasov 2004). However, this definition of harmonic coordinates then becomes theory-dependent (Will 2018b). In order to alleviate the ambiguity, Kopeikin and Vlasov (2004); Kopeikin (2019) used the name “Nutku gauge”—or “conformal harmonic gauge”—when the harmonic conditions are imposed on the metric in the Einstein frame rather than in the Jordan frame. Will (2018b), it corresponds to the “generalized harmonic gauge”.

### 4.1.1 Equations of motion, Lense–Thirring and Shapiro delay

In harmonic coordinates, the EIHDL equation of motion in the barycentric frame [Eq. (90)] reads as follows<fnref n="25" />

<footnote n="25">²⁵ However, the transformations between the barycentric and geocentric reference frames slightly depend on the specific type of harmonic coordinates being used Kopeikin and Vlasov (2004).</footnote>
$$ 
\begin{aligned}
g_{00} &= -1 + \frac{2w}{c^2} - \frac{2\beta w^2}{c^4} + \mathcal{O}(c^{-5}) \\
g_{0i} &= -\frac{(4\gamma + 3)w_i}{2c^3} + \mathcal{O}(c^{-4}) \\
g_{ij} &= \delta_{ij} \left(1 + \frac{2\gamma w}{c^2}\right) + \mathcal{O}(c^{-3}),
\end{aligned} 
\quad (95) $$

In general relativity, the PPN parameters $\beta$ and $\gamma$ equal to one, but not for alternative theories in general.

The form of the metric given in Eq. (95) is not sufficient in order to completely fix the coordinate system, and one has to impose an additional constraint (or gauge) in order to specify the metric field equations—see, e.g., Klioner and Soffel (2000). Since the IAU recommends the coordinates (or gauge) to be harmonic for general relativity (Soffel 2003), it seems relevant to consider harmonic coordinates for alternative theories as well, such that (at least) the coordinate systems of the alternative theory reduce to the IAU coordinate systems in the limit for which both theories cannot be distinguished at the level of accuracy of the observations.

### Ambiguity of the definition of harmonic coordinates

The very definition of what one may call harmonic coordinates can be ambiguous for alternative theories of gravity. For instance, for scalar-tensor theories, one can impose the usual post-Newtonian harmonic conditions—$g^{\alpha\beta}\Gamma^0_{\alpha\beta} = 0$—either in the Jordan frame, or in the conformal Einstein frame (Damour and Esposito-Farese 1992; Kopeikin and Vlasov 2004; Minazzoli and Chauvineau 2011; Will 2018b; Kopeikin 2019). But both choices reduce to the general relativity harmonic coordinate in the limit for which the difference between general relativity and the scalar-tensor theory becomes negligible (e.g. for the Brans–Dicke parameter that is such that $\omega \to \infty$). Imposing the harmonic conditions in the Einstein frame has been shown to lead to more convenient properties (Damour and Esposito-Farese 1992; Kopeikin and Vlasov 2004; Minazzoli and Chauvineau 2011; Kopeikin 2019), notably with respect to the full theory of reference frames (Kopeikin and Vlasov 2004). However, this definition of harmonic coordinates then becomes theory-dependent (Will 2018b). In order to alleviate the ambiguity, Kopeikin and Vlasov (2004); Kopeikin (2019) used the name “Nutku gauge”—or “conformal harmonic gauge”—when the harmonic conditions are imposed on the metric in the Einstein frame rather than in the Jordan frame. Will (2018b), it corresponds to the “generalized harmonic gauge”.

### 4.1.1 Equations of motion, Lense–Thirring and Shapiro delay

In harmonic coordinates, the EIHDL equation of motion in the barycentric frame [Eq. (90)] reads as follows<fnref n="25" />

<footnote n="25">²⁵ However, the transformations between the barycentric and geocentric reference frames slightly depend on the specific type of harmonic coordinates being used Kopeikin and Vlasov (2004).</footnote>
$$
\begin{equation}
\begin{aligned}
\boldsymbol{a}_T^{\text{PPN}} = & -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3} \boldsymbol{r}_{AT} \\
& -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3 c^2} \boldsymbol{r}_{AT} \left\{ \gamma v_T^2 + (\gamma + 1)v_A^2 - 2(1 + \gamma)\boldsymbol{v}_A \cdot \boldsymbol{v}_T - \frac{3}{2} \left( \frac{\boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A}{r_{AT}} \right)^2 \right. \\
& \left. - \frac{1}{2} \boldsymbol{r}_{AT} \cdot \boldsymbol{a}_A - 2(\gamma + \beta) \sum_{B \neq T} \frac{\mu_B}{r_{TB}} - (2\beta - 1) \sum_{B \neq A} \frac{\mu_B}{r_{AB}} \right\} \\
& + \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}^3} \left[ 2(1 + \gamma) \boldsymbol{r}_{AT} \cdot \boldsymbol{v}_T - (1 + 2\gamma) \boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A \right] (\boldsymbol{v}_T - \boldsymbol{v}_A) \\
& + \frac{3 + 4\gamma}{2} \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}} \boldsymbol{a}_A,
\end{aligned}
\tag{96}
\end{equation}
$$

with $\gamma$ and $\beta$, the Parametrized Post Newtonian (PPN) parameters. They are equal to 1 in general relativity, such that one recovers the usual EIHDL Eq. (90) in that case.

Additionally, directly from the metric Eq. (95), one can deduce that the Shapiro delay becomes

$$
c(t_r - t_e)^{\text{PPN}} = R + \sum_A (\gamma + 1) \frac{\mu_A}{c^2} \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}, \quad (97)
$$

with the same notations as Eq. (94).

The Lense–Thirring acceleration is also modified and reads:

$$
\boldsymbol{a}_{LT}^{\text{PPN}} = (1 + \gamma) \frac{GS}{c^2 r^3} \left[ \frac{3\boldsymbol{k} \cdot \boldsymbol{r}}{r^2} (\boldsymbol{r} \wedge \boldsymbol{v}) - (\boldsymbol{k} \wedge \boldsymbol{v}) \right] \quad (98)
$$

with the notation of Eq. (91).

It is also important to stress that only the PPN parameter $\gamma$ appears in the Shapiro delay formula that is used for the computation of the propagation time of light and doppler shift (see Sect. 3.4). As a consequence, the impact of this parameter can be more easily disentangled from the effect of other parameters, in particular with conjunction events for which one directly probes the differential evolution of Eq. (97) with respect to evolution of the impact parameter during the conjunction (Bertotti et al. 2003b). This disentangling of the $\gamma$ parameter is also enhanced by the LT contribution, which depends on $\gamma$.

Since the first radar echos obtained from the Mercury surface (Shapiro 1964), it
has been possible to evaluate the departure from unity of the PPN parameters in the
context of the planetary orbit computation. Such estimations are done during the
construction of the planetary ephemerides when the initial conditions of the planetary
orbits are fitted to observations with a least square procedure (Ash et al. 1967;
Standish et al. 1976), together with other parameters such asteroid masses.
$$ 
\begin{aligned}
\boldsymbol{a}_{T}^{\text{PPN}} = & -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3} \boldsymbol{r}_{AT} \\
& -\sum_{A \neq T} \frac{\mu_A}{r_{AT}^3 c^2} \boldsymbol{r}_{AT} \left\{ \gamma v_T^2 + (\gamma + 1)v_A^2 - 2(1 + \gamma)\boldsymbol{v}_A \cdot \boldsymbol{v}_T - \frac{3}{2} \left( \frac{\boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A}{r_{AT}} \right)^2 \right. \\
& \left. - \frac{1}{2} \boldsymbol{r}_{AT} \cdot \boldsymbol{a}_A - 2(\gamma + \beta) \sum_{B \neq T} \frac{\mu_B}{r_{TB}} - (2\beta - 1) \sum_{B \neq A} \frac{\mu_B}{r_{AB}} \right\} \\
& + \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}^3} \left[ 2(1 + \gamma) \boldsymbol{r}_{AT} \cdot \boldsymbol{v}_T - (1 + 2\gamma) \boldsymbol{r}_{AT} \cdot \boldsymbol{v}_A \right] (\boldsymbol{v}_T - \boldsymbol{v}_A) \\
& + \frac{3 + 4\gamma}{2} \sum_{A \neq T} \frac{\mu_A}{c^2 r_{AT}} \boldsymbol{a}_A,
\end{aligned}
\quad (96) $$

with $\gamma$ and $\beta$, the Parametrized Post Newtonian (PPN) parameters. They are equal to 1 in general relativity, such that one recovers the usual EIHDL Eq. (90) in that case.

Additionally, directly from the metric Eq. (95), one can deduce that the Shapiro delay becomes

$$ c(t_r - t_e)^{\text{PPN}} = R + \sum_A (\gamma + 1) \frac{\mu_A}{c^2} \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}, \quad (97) $$

with the same notations as Eq. (94).

The Lense–Thirring acceleration is also modified and reads:

$$ \boldsymbol{a}_{LT}^{\text{PPN}} = (1 + \gamma) \frac{GS}{c^2 r^3} \left[ \frac{3\boldsymbol{k} \cdot \boldsymbol{r}}{r^2} (\boldsymbol{r} \wedge \boldsymbol{v}) - (\boldsymbol{k} \wedge \boldsymbol{v}) \right] \quad (98) $$

with the notation of Eq. (91).

It is also important to stress that only the PPN parameter $\gamma$ appears in the Shapiro delay formula that is used for the computation of the propagation time of light and doppler shift (see Sect. 3.4). As a consequence, the impact of this parameter can be more easily disentangled from the effect of other parameters, in particular with conjunction events for which one directly probes the differential evolution of Eq. (97) with respect to evolution of the impact parameter during the conjunction (Bertotti et al. 2003b). This disentangling of the $\gamma$ parameter is also enhanced by the LT contribution, which depends on $\gamma$.

Since the first radar echos obtained from the Mercury surface (Shapiro 1964), it has been possible to evaluate the departure from unity of the PPN parameters in the context of the planetary orbit computation. Such estimations are done during the construction of the planetary ephemerides when the initial conditions of the planetary orbits are fitted to observations with a least square procedure (Ash et al. 1967; Standish et al. 1976), together with other parameters such asteroid masses.
## 4.1.2 Time-scales

As mentioned in Sect. 3.2.2, planetary ephemerides produce the difference between the TT and the TDB. In the PPN framework, Eq. (86) has to be modified as follows (Manche 2011)

$$ 
\begin{aligned}
a &= -\frac{1}{2}v_T^2 - \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \\
b &= -\frac{1}{8}v_T^4 + \left(\beta - \frac{1}{2}\right) \left[ \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \right]^2 + \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \left\{ \begin{aligned} &2(1 + \gamma)v_T \cdot \mathbf{v}_A - (\gamma + \frac{1}{2})v_T^2 \\ &-(1 + \gamma)v_A^2 + \frac{1}{2}\mathbf{a}_A \cdot \mathbf{r}_{AT} + \frac{1}{2}\left(\frac{\mathbf{v}_A \cdot \mathbf{r}_{AT}}{r_{AT}}\right)^2 \end{aligned} \right\} \\ 
&+ (2\beta - 1) \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \left( \sum_{B \neq A} \frac{\mu_B}{r_{BA}} \right).
\end{aligned} 
\quad (99) $$

### Comments on reference frame theory in alternative theories

The difficult subject of reference frame systems in alternative theories has not been investigated systematically. In fact, it has been explored explicitly only for a class of massless scalar-tensor theories “à la Brans–Dicke” (Kopeikin and Vlasov 2004), after some previous effort to characterize it in a general standard post-Newtonian approach (Kloner and Soffel 2000). Fortunately in practice, the differences between coordinate systems in general relativity versus PPN framework are numerically negligible—notably given the already tight constraints that one already has on the PPN parameters.

## 4.1.3 Constraints and correlations

An important discussion related to the PPN parameter determinations using least square methods is the correlations between $\beta$, $\gamma$ and other contributions to the planetary accelerations. This problem has been pointed out a long time ago (Anderson et al. 1978) but it is still vivid today. In particular, if we consider the analytical expression of the advance of the perihelia for planetary orbits $\Delta\dot{\omega}_{PLA}$, as given in Will (2014), the two main contributions of this equation are the PPN term $(2 + 2\gamma - \beta)$ and a term giving the impact of the Sun oblateness $J_2^\odot$. But we can also add the main belt asteroid contribution and the LT effect defined in Sect. 3.3.3, such as per orbit, one gets (Bertotti et al. 2003a)

$$ 
\begin{aligned}
\Delta\dot{\omega}_{PLA} = \mu_{\odot}(t) &\times \left( \frac{2(2\gamma - \beta + 2)}{3a(1 - e^2)c^2} + \Delta\dot{\omega}_{LT}(S, \gamma) \right) \\ 
&+ J_2^\odot \frac{R_{\odot}^2}{a^2(1 - e^2)^2} + \Delta\dot{\omega}_{AST}(\mu_{AST}, a_{AST}, e_{AST})
\end{aligned} 
\quad (100) $$

where $a_{AST}$, $e_{AST}$, $\mu_{AST}$ are respectively the asteroid orbit semi-major axis and eccentricity and its gravitational mass perturbing the planet orbits (main-belt and
### 4.1.2 Time-scales

As mentioned in Sect. 3.2.2, planetary ephemerides produce the difference between the TT and the TDB. In the PPN framework, Eq. (86) has to be modified as follows (Manche 2011)

$$ 
\begin{aligned}
a &= -\frac{1}{2}v_T^2 - \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \\
b &= -\frac{1}{8}v_T^4 + \left(\beta - \frac{1}{2}\right) \left[ \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \right]^2 + \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \left\{ \begin{aligned} &2(1 + \gamma)v_T \cdot v_A - (\gamma + \frac{1}{2})v_T^2 \\ &- (1 + \gamma)v_A^2 + \frac{1}{2}a_A \cdot r_{AT} + \frac{1}{2} \left( \frac{v_A \cdot r_{AT}}{r_{AT}} \right)^2 \end{aligned} \right\} \\ 
&\quad + (2\beta - 1) \sum_{A \neq T} \frac{\mu_A}{r_{AT}} \left( \sum_{B \neq A} \frac{\mu_B}{r_{BA}} \right).
\end{aligned} 
\quad (99) $$

*Comments on reference frame theory in alternative theories*

The difficult subject of reference frame systems in alternative theories has not been investigated systematically. In fact, it has been explored explicitly only for a class of massless scalar-tensor theories “à la Brans–Dicke” (Kopeikin and Vlasov 2004), after some previous effort to characterize it in a general standard post-Newtonian approach (Kloner and Soffel 2000). Fortunately in practice, the differences between coordinate systems in general relativity versus PPN framework are numerically negligible—notably given the already tight constraints that one already has on the PPN parameters.

### 4.1.3 Constraints and correlations

An important discussion related to the PPN parameter determinations using least square methods is the correlations between $\beta$, $\gamma$ and other contributions to the planetary accelerations. This problem has been pointed out a long time ago (Anderson et al. 1978) but it is still vivid today. In particular, if we consider the analytical expression of the advance of the perihelia for planetary orbits $\Delta\dot{\omega}_{PLA}$, as given in Will (2014), the two main contributions of this equation are the PPN term $(2 + 2\gamma - \beta)$ and a term giving the impact of the Sun oblateness $J_2^\odot$. But we can also add the main belt asteroid contribution and the LT effect defined in Sect. 3.3.3, such as per orbit, one gets (Bertotti et al. 2003a)

$$ 
\begin{aligned}
\Delta\dot{\omega}_{PLA} = \mu_{\odot}(t) \times & \left( \frac{2(2\gamma - \beta + 2)}{3a(1 - e^2)c^2} + \Delta\dot{\omega}_{LT}(S, \gamma) \right) \\ 
& + J_2^\odot \frac{R_{\odot}^2}{a^2(1 - e^2)^2} + \Delta\dot{\omega}_{AST}(\mu_{AST}, a_{AST}, e_{AST})
\end{aligned} 
\quad (100) $$

where $a_{AST}$, $e_{AST}$, $\mu_{AST}$ are respectively the asteroid orbit semi-major axis and eccentricity and its gravitational mass perturbing the planet orbits (main-belt and
trans-neptunian objects), a, e are the semi-major axis and eccentricity of the planet
orbit and finally, c and µ⊙ are the speed of light and the gravitational parameter of
the Sun. This latest quantity can be time-varying as discussed in Sect. 4.3.
Δω$_{LT}$ (S, γ) is the Lense–Thirring effect presented in Sect. 3.3.3. A discussion about
other contributing terms such as PN cross-terms, de Sitter precession and 2PN
crossterms has been proposed by Will (2018). The conclusion is that these 3 addi-
tional contributions could have some importances in the coming years with the Bepi-
Colombo mission. In any case, they are implicitly taken into account when
numerically integrating the EIHDL equation of motion.

Note that Eq. (100) is an analytical approximation of the perihelia advance rate that is useful in order to get a sense of some of the correlations that may arise between different quantities. The actual value would have to be inferred from the numerical integration of the equations of motion in the solar system with the total acceleration, **a**<sub>T</sub><sup>total:</sup>:

$$
\mathbf{a}_T^{\text{total}} = \mathbf{a}_T(\mu_{SS}, \gamma, \beta, IC_{SS}) + \mathbf{a}_{J_2^\odot}(\mu_\odot, J_2^\odot) + \mathbf{a}_{LT}(S_\odot, \gamma), \quad (101)
$$

where **a**<sub>T</sub>, **a**<sub>J₂</sub><sup>☉</sup>, **a**<sub>LT</sub> are the point-mass interaction acceleration given in Eq. (96), the
acceleration induced by the oblateness of the Sun Eq. (92) and the Lense-Thirring
acceleration Eq. (98) respectively, *IC*<sub>SS</sub> stands for the initial conditions of the
numerical integration for the solar system bodies, including asteroids.

It is rather understandable from Eq. (100) that there can be some level of
correlation between, e.g., the semi-major axis (which follows from the initial
conditions in the numerical integration) and the post-Newtonian parameters γ and β.
One can also see that for disentangling the different contributions, it is more efficient
to have several orbits with significantly different semi-major axis a and eccentricities
e as the PPN acceleration depends on a when the Sun oblateness acceleration
depends on a² [see Eqs. (100) and (93)]. It is also interesting to compare the Sun J₂⊙
determinations obtained in general relativity—equivalent to consider β and γ fixed to
one—presented on Fig. 6 with those obtained in adding the PPN parameters β
and γ into the global planetary fit as on Fig. 9. It appears first that, besides the
differences between ephemerides, the general relativity estimations (with β and γ
fixed to unity) have globally reduced error bars (i.e. less uncertainty induced by less
correlation), as one can see, for example, with INPOP20a (2.218 ± 0.03 in general
relativity versus 2.165 ± 0.12 in PPN), or the partial fit obtained by Genova et al.
(2018) (2.271 ± 0.003 in general relativity versus 2.246 ± 0.02 in PPN). Secondly,
biases affect also the J₂⊙ determinations obtained with fitted PPN parameters as these
latest can be significantly different from the J₂⊙ values obtained in general relativity.
These bias are due to the correlations between the PPN parameters and J₂⊙. In the
case of the Earth–Moon system, the relation is even more complex as it also includes
the tidal contributions depending on the internal structures of the two bodies.

In order to overcome the problem of the correlations between the parameters,
other methods of inversion such as random walk exploration or Genetic Algorithm
have been proposed (see Sect. 3.5). They usually give more conservative constraints
than the direct least square method as it is visible on Table 5 and Fig. 8. On these
trans-neptunian objects), a, e are the semi-major axis and eccentricity of the planet
orbit and finally, c and µ⊙ are the speed of light and the gravitational parameter of
the Sun. This latest quantity can be time-varying as discussed in Sect. 4.3.
Δω$_{LT}$ (S, γ) is the Lense–Thirring effect presented in Sect. 3.3.3. A discussion about
other contributing terms such as PN cross-terms, de Sitter precession and 2PN
crossterms has been proposed by Will (2018). The conclusion is that these 3 addi-
tional contributions could have some importances in the coming years with the Bepi-
Colombo mission. In any case, they are implicitly taken into account when
numerically integrating the EIHDL equation of motion.

Note that Eq. (100) is an analytical approximation of the perihelia advance rate that is useful in order to get a sense of some of the correlations that may arise between different quantities. The actual value would have to be inferred from the numerical integration of the equations of motion in the solar system with the total acceleration, **a**<sub>T</sub><sup>total:</sup>

$$
\mathbf{a}_T^{\text{total}} = \mathbf{a}_T(\mu_{SS}, \gamma, \beta, IC_{SS}) + \mathbf{a}_{J_2^\odot}(\mu_\odot, J_2^\odot) + \mathbf{a}_{LT}(S_\odot, \gamma), \quad (101)
$$

where **a**<sub>T</sub>, **a**<sub>J₂</sub><sup>☉</sup>, **a**<sub>LT</sub> are the point-mass interaction acceleration given in Eq. (96), the
acceleration induced by the oblateness of the Sun Eq. (92) and the Lense-Thirring
acceleration Eq. (98) respectively, *IC*<sub>SS</sub> stands for the initial conditions of the
numerical integration for the solar system bodies, including asteroids.

It is rather understandable from Eq. (100) that there can be some level of
correlation between, e.g., the semi-major axis (which follows from the initial
conditions in the numerical integration) and the post-Newtonian parameters γ and β.
One can also see that for disentangling the different contributions, it is more efficient
to have several orbits with significantly different semi-major axis a and eccentricities
e as the PPN acceleration depends on a when the Sun oblateness acceleration
depends on a² [see Eqs. (100) and (93)]. It is also interesting to compare the Sun J₂⊙
determinations obtained in general relativity—equivalent to consider β and γ fixed to
one—presented on Fig. 6 with those obtained in adding the PPN parameters β
and γ into the global planetary fit as on Fig. 9. It appears first that, besides the
differences between ephemerides, the general relativity estimations (with β and γ
fixed to unity) have globally reduced error bars (i.e. less uncertainty induced by less
correlation), as one can see, for example, with INPOP20a (2.218 ± 0.03 in general
relativity versus 2.165 ± 0.12 in PPN), or the partial fit obtained by Genova et al.
(2018) (2.271 ± 0.003 in general relativity versus 2.246 ± 0.02 in PPN). Secondly,
biases affect also the J₂⊙ determinations obtained with fitted PPN parameters as these
latest can be significantly different from the J₂⊙ values obtained in general relativity.
These bias are due to the correlations between the PPN parameters and J₂⊙. In the
case of the Earth–Moon system, the relation is even more complex as it also includes
the tidal contributions depending on the internal structures of the two bodies.

In order to overcome the problem of the correlations between the parameters,
other methods of inversion such as random walk exploration or Genetic Algorithm
have been proposed (see Sect. 3.5). They usually give more conservative constraints
than the direct least square method as it is visible on Table 5 and Fig. 8. On these
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Constraints on PPN parameters $\beta$ and $\gamma$ obtained from planetary ephemerides as presented in Table 5. The shadded area is the zone common to all the constraints. The lightgrey lines give the constraint on $\gamma$ extracted from the Cassini experiment (Bertotti et al. 2003b). Three categories are considered: the values obtained from a full global fit, the one obtained using random walk exploration methods (RW) and finally, the one constrained either in fixing one of the two parameters [usually $\gamma$ like in Konopliv et al. (2011); Park et al. (2017)] either in considering only one planetary orbit (like the Mercury orbit for Genova et al. (2018). These categories were labelled according to the Column 3 of Table 5 and the annotated text refers to Column 2 of Table 5</figcaption></figure>

table and figure, one can find the estimations of PPN parameters $\beta$ and $\gamma$ determinations based on planetary ephemerides. Three families of determinations are proposed. The first noted *FF* gathers least square fitted evaluations of PPN parameters, $J_2^\circ$ and secular variation of the gravitational parameter of the Sun, $\mu_\odot$. As discussed previously on Sects. 3.3.4 and 3.5, these parameters are strongly correlated and thus are affected by biases (see Figs. 6, 9) and underestimated uncertainties. The second type of determinations given on Table 5 with the label *PF* are the one deduced with one (at least) of these parameters fixed in the planetary ephemeris adjustment. These estimations are less correlated than the *FF* ones and give more reliable assessments but are limited to only one value (usually either $\gamma$ or $J_2^\circ$ are fixed and $\beta$ is fitted). Finally, the latest category is obtained by Genetic Algorithm (Fienga et al. 2015) or random walk exploration (Fienga et al. 2021a). They investigate the space of parameters more largely than the two former methods, leading to usually larger intervals for possible general relativity violations. On Table 5, values obtained in the context of one single planet analysis (Mercury) with only one single spacecraft (MESSENGER) data analysis (Genova et al. 2018) are also given for comparisons. The obtained results show smaller uncertainties than the planetary determinations but might be affected by possible biases. As it was shown in the Eq. (100), the use of
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Constraints on PPN parameters $\beta$ and $\gamma$ obtained from planetary ephemerides as presented in Table 5. The shadded area is the zone common to all the constraints. The lightgrey lines give the constraint on $\gamma$ extracted from the Cassini experiment (Bertotti et al. 2003b). Three categories are considered: the values obtained from a full global fit, the one obtained using random walk exploration methods (RW) and finally, the one constrained either in fixing one of the two parameters [usually $\gamma$ like in Konopliv et al. (2011); Park et al. (2017)] either in considering only one planetary orbit (like the Mercury orbit for Genova et al. (2018). These categories were labelled according to the Column 3 of Table 5 and the annotated text refers to Column 2 of Table 5</figcaption></figure>

table and figure, one can find the estimations of PPN parameters $\beta$ and $\gamma$ determinations based on planetary ephemerides. Three families of determinations are proposed. The first noted *FF* gathers least square fitted evaluations of PPN parameters, $J_2^\circ$ and secular variation of the gravitational parameter of the Sun, $\mu_\odot$. As discussed previously on Sects. 3.3.4 and 3.5, these parameters are strongly correlated and thus are affected by biases (see Figs. 6, 9) and underestimated uncertainties. The second type of determinations given on Table 5 with the label *PF* are the one deduced with one (at least) of these parameters fixed in the planetary ephemeris adjustment. These estimations are less correlated than the *FF* ones and give more reliable assessments but are limited to only one value (usually either $\gamma$ or $J_2^\circ$ are fixed and $\beta$ is fitted). Finally, the latest category is obtained by Genetic Algorithm (Fienga et al. 2015) or random walk exploration (Fienga et al. 2021a). They investigate the space of parameters more largely than the two former methods, leading to usually larger intervals for possible general relativity violations. On Table 5, values obtained in the context of one single planet analysis (Mercury) with only one single spacecraft (MESSENGER) data analysis (Genova et al. 2018) are also given for comparisons. The obtained results show smaller uncertainties than the planetary determinations but might be affected by possible biases. As it was shown in the Eq. (100), the use of

<figure><img src="image_12.png" /></figure>
<table><caption>**Table 5** PPN parameters $\beta$ and $\gamma$ obtained with planetary ephemerides. Columns 1, 2 and 3 give respectively the reference, the ephemeris involved and method used (MSG). As explained in the text (Sect. 4.1), we consider a classification based on 3 types of adjustment methods: the direct least square fit combining $\beta$, $\gamma$, all the ephemeris parameters and the oblateness of the Sun, $J_2^\odot$ (FF), the partial fit where some parameters are fixed such as $\gamma$ or $J_2^\odot$ and the other are fitted (PF), the random walk exploration (RW) where $\beta$ and $\gamma$ are randomly sampled when the rest of the parameters (including $J_2^\odot$) are fitted. The results obtained solely from considering Mercury orbit (MSG) are provided for informational purposes</caption><thead><tr><th>References</th><th>PE</th><th>Method</th><th>β − 1 × 10<sup>5</sup></th><th>γ − 1 × 10<sup>5</sup></th><th>J<sub>2</sub><sup>☉</sup> × 10<sup>7</sup></th></tr></thead><tbody><tr><td>Pitjeva (2001)</td><td>EPM2000</td><td>FF</td><td>40 ± 20*</td><td>10 ± 10*</td><td>2.43 ± 0.67*</td></tr><tr><td>Standish (2001)</td><td>DE405</td><td>FF</td><td>10 ± 10*</td><td>40 ± 10*</td><td>2.46 ± 0.68*</td></tr><tr><td>Pitjeva (2005b)</td><td>EPM2004</td><td>FF</td><td>0 ± 10*</td><td>−10 ± 20*</td><td>1.9 ± 0.3</td></tr><tr><td>Fienga et al. (2009)</td><td>INPOP08a</td><td>PF</td><td>7.5 ± 12.5</td><td>0.0</td><td>fixed to 1.82</td></tr><tr><td>Fienga et al. (2011a)</td><td>INPOP10a</td><td>PF</td><td>−4.1 ± 7.8</td><td>−6.2 ± 8.1</td><td>fixed to 2.40</td></tr><tr><td>Konopliv et al. (2011)</td><td>DE421</td><td>PF</td><td>4 ± 24</td><td>fixed to 2.1</td><td>fixed to 1.8</td></tr><tr><td>Pitjeva and Pitjev (2013)</td><td>EPM2011</td><td>FF</td><td>−2 ± 3</td><td>4 ± 6</td><td>2.0 ± 0.20</td></tr><tr><td>Verma et al. 2014</td><td>INPOP13c</td><td>FF</td><td>0.2 ± 2.5</td><td>−0.3 ± 2.5</td><td>2.40 ± 0.20</td></tr><tr><td>Fienga et al. (2015)</td><td>INPOP15a</td><td>FF</td><td>−6.7 ± 6.9</td><td>−0.8 ± 5.7</td><td>2.27 ± 0.25</td></tr><tr><td>Fienga et al. (2015)</td><td>INPOP15a</td><td>RW</td><td>0.00 ± 6.90</td><td>−1.55 ± 5.01</td><td>2.22 ± 0.13</td></tr><tr><td>Park et al. (2017)</td><td>DE440</td><td>PF</td><td>−2.6 ± 3.9*</td><td>fixed to 2.1</td><td>2.25 ± 0.09*</td></tr><tr><td>Genova et al. (2018)</td><td>DE438</td><td>MSG</td><td>−1.625 ± 1.8</td><td>fixed to 2.1</td><td>2.246 ± 0.02</td></tr><tr><td>Fienga et al. (2021a)</td><td>INPOP20a</td><td>PF</td><td>1.9 ± 6.28</td><td>2.64 ± 3.44</td><td>2.165 ± 0.12</td></tr><tr><td>Fienga et al. (2021a)</td><td>INPOP20a</td><td>RW</td><td>−1.12 ± 7.16</td><td>−1.69 ± 7.49</td><td>2.206 ± 0.03</td></tr></tbody></table>

multiple orbits favors the disentangling of the different planetary contributions. Moreover, the value of $\beta$ is obtained by introducing the Nordvedt relation Eq. (126) in the equation of motion and in the fitting procedure which supposes to consider only metric theories and increases the correlations between the fitted parameters (see Sect. 4.2).

Finally, as one can see on Fig. 8, all the PPN estimations obtained with the three methods are consistent with general relativity, with no significant departure from unity for $\beta$ and $\gamma$. By considering the common overlaps of all the recently published intervals obtained with planetary ephemerides, the two following constraints are obtained:

$$ (\beta - 1) = (-0.45 \pm 1.75) \times 10^{-5} \\ (\gamma - 1) = (0.55 \pm 1.35) \times 10^{-5} \qquad (102) $$

These values represent an improvement of about 3 orders of magnitude for $\gamma$ and $\beta$ in comparison to the historic first determination of Anderson et al. (1978). We can expect the same type of improvements in the coming years, especially with the outcomes of the BepiColombo mission (Imperi et al. 2018; De Marchi and Cascioli 2020).
<table><caption>**Table 5** PPN parameters $\beta$ and $\gamma$ obtained with planetary ephemerides. Columns 1, 2 and 3 give respectively the reference, the ephemeris involved and method used (MSG). As explained in the text (Sect. 4.1), we consider a classification based on 3 types of adjustment methods: the direct least square fit combining $\beta$, $\gamma$, all the ephemeris parameters and the oblateness of the Sun, $J_2^\odot$ (FF), the partial fit where some parameters are fixed such as $\gamma$ or $J_2^\odot$ and the other are fitted (PF), the random walk exploration (RW) where $\beta$ and $\gamma$ are randomly sampled when the rest of the parameters (including $J_2^\odot$) are fitted. The results obtained solely from considering Mercury orbit (MSG) are provided for informational purposes</caption><thead><tr><th>References</th><th>PE</th><th>Method</th><th>β − 1 × 10<sup>5</sup></th><th>γ − 1 × 10<sup>5</sup></th><th>J<sub>2</sub><sup>☉</sup> × 10<sup>7</sup></th></tr></thead><tbody><tr><td>Pitjeva (2001)</td><td>EPM2000</td><td>FF</td><td>40 ± 20*</td><td>10 ± 10*</td><td>2.43 ± 0.67*</td></tr><tr><td>Standish (2001)</td><td>DE405</td><td>FF</td><td>10 ± 10 *</td><td>40 ± 10*</td><td>2.46 ± 0.68*</td></tr><tr><td>Pitjeva (2005b)</td><td>EPM2004</td><td>FF</td><td>0 ± 10 *</td><td>−10 ± 20 *</td><td>1.9 ± 0.3</td></tr><tr><td>Fienga et al. (2009)</td><td>INPOP08a</td><td>PF</td><td>7.5 ± 12.5</td><td>0.0</td><td>fixed to 1.82</td></tr><tr><td>Fienga et al. (2011a)</td><td>INPOP10a</td><td>PF</td><td>−4.1 ± 7.8</td><td>−6.2 ± 8.1</td><td>fixed to 2.40</td></tr><tr><td>Konopliv et al. (2011)</td><td>DE421</td><td>PF</td><td>4 ± 24</td><td>fixed to 2.1</td><td>fixed to 1.8</td></tr><tr><td>Pitjeva and Pitjev (2013)</td><td>EPM2011</td><td>FF</td><td>−2 ± 3</td><td>4 ± 6</td><td>2.0 ± 0.20</td></tr><tr><td>Verma et al. 2014</td><td>INPOP13c</td><td>FF</td><td>0.2 ± 2.5</td><td>−0.3 ± 2.5</td><td>2.40 ± 0.20</td></tr><tr><td>Fienga et al. (2015)</td><td>INPOP15a</td><td>FF</td><td>−6.7 ± 6.9</td><td>−0.8 ± 5.7</td><td>2.27 ± 0.25</td></tr><tr><td>Fienga et al. (2015)</td><td>INPOP15a</td><td>RW</td><td>0.00 ± 6.90</td><td>−1.55 ± 5.01</td><td>2.22 ± 0.13</td></tr><tr><td>Park et al. (2017)</td><td>DE440</td><td>PF</td><td>−2.6 ± 3.9*</td><td>fixed to 2.1</td><td>2.25 ± 0.09*</td></tr><tr><td>Genova et al. (2018)</td><td>DE438</td><td>MSG</td><td>−1.625 ± 1.8</td><td>fixed to 2.1</td><td>2.246 ± 0.02</td></tr><tr><td>Fienga et al. (2021a)</td><td>INPOP20a</td><td>PF</td><td>1.9 ± 6.28</td><td>2.64 ± 3.44</td><td>2.165 ± 0.12</td></tr><tr><td>Fienga et al. (2021a)</td><td>INPOP20a</td><td>RW</td><td>−1.12 ± 7.16</td><td>−1.69 ± 7.49</td><td>2.206 ± 0.03</td></tr></tbody></table>

multiple orbits favors the disentangling of the different planetary contributions. Moreover, the value of $\beta$ is obtained by introducing the Nordvedt relation Eq. (126) in the equation of motion and in the fitting procedure which supposes to consider only metric theories and increases the correlations between the fitted parameters (see Sect. 4.2).

Finally, as one can see on Fig. 8, all the PPN estimations obtained with the three methods are consistent with general relativity, with no significant departure from unity for $\beta$ and $\gamma$. By considering the common overlaps of all the recently published intervals obtained with planetary ephemerides, the two following constraints are obtained:

$$ (\beta - 1) = (-0.45 \pm 1.75) \times 10^{-5} \\ (\gamma - 1) = (0.55 \pm 1.35) \times 10^{-5} \qquad (102) $$

These values represent an improvement of about 3 orders of magnitude for $\gamma$ and $\beta$ in comparison to the historic first determination of Anderson et al. (1978). We can expect the same type of improvements in the coming years, especially with the outcomes of the BepiColombo mission (Imperi et al. 2018; De Marchi and Cascioli 2020).
<figure><img src="image_11.png" /><figcaption>Fig. 9 Values of the Sun oblateness $J_2^\circ$ obtained when considering PPN parameters $\beta$ and $\gamma$ as varying parameters. The dotted lines give the limits of the less constraining helioseismologic value from Pijpers (1998). F stands for Full fit, RW for random walk exploration method and P for partially fitted ephemerides (some parameters are fixed). N (respectively Y) indicates that the value has been obtained without (respectively with) LT. In y-axis, are given the references: G18 is Genova et al. (2018), I20a (Fienga et al. 2021a), D440 (Park et al. 2017), I13c (Fienga et al. 2014), E04 and E11 in Pitjeva and Pitjev (2013). E00 is Pitjeva (2001), D405 (Standish 2001) and I15a (Fienga et al. 2015). It is important to note that in Genova et al. (2018), only the orbit of Mercury is considered</figcaption></figure>

### 4.1.4 The extended PPN formalism

#### SPN gauge versus harmonic gauge

Usually, the extended standard post-Newtonian formalism assumes a gauge (or coordinate system) that is quite different from the one recommended by the IAU that is called the standard post-Newtonian (SPN) gauge. In general relativity already, while the harmonic condition<fnref n="26" /> at the post-Newtonian level reads

$$g^{\alpha\beta} \Gamma^{0}_{\alpha\beta} = \mathcal{O}(c^{-5}), \qquad (103)$$

the SPN conditions read (Damour et al. 1991)

$$\hat{\partial}_j g_{0j} - \frac{1}{2} \hat{\partial}_0 g_{jj} = \mathcal{O}(c^{-5}), \qquad (104)$$

<footnote n="26">²⁶ On top of the strong isotropy condition imposed on the metric, which fixes the space coordinate (Damour et al. 1991; Soffel 2003)—see Sect. 2.1.7.</footnote>
<figure><img src="image_13.png" /><figcaption>Fig. 9 Values of the Sun oblateness $J_2^\odot$ obtained when considering PPN parameters $\beta$ and $\gamma$ as varying parameters. The dotted lines give the limits of the less constraining helioseismologic value from Pijpers (1998). F stands for Full fit, RW for random walk exploration method and P for partially fitted ephemerides (some parameters are fixed). N (respectively Y) indicates that the value has been obtained without (respectively with) LT. In y-axis, are given the references: G18 is Genova et al. (2018), I20a (Fienga et al. 2021a), D440 (Park et al. 2017), I13c (Fienga et al. 2014), E04 and E11 in Pitjeva and Pitjev (2013). E00 is Pitjeva (2001), D405 (Standish 2001) and I15a (Fienga et al. 2015). It is important to note that in Genova et al. (2018), only the orbit of Mercury is considered</figcaption></figure>

### 4.1.4 The extended PPN formalism

#### SPN gauge versus harmonic gauge

Usually, the extended standard post-Newtonian formalism assumes a gauge (or coordinate system) that is quite different from the one recommended by the IAU that is called the standard post-Newtonian (SPN) gauge. In general relativity already, while the harmonic condition<fnref n="26" /> at the post-Newtonian level reads

$$g^{\alpha\beta} \Gamma^{0}_{\alpha\beta} = \mathcal{O}(c^{-5}), \qquad (103)$$

the SPN conditions read (Damour et al. 1991)

$$\hat{\partial}_j g_{0j} - \frac{1}{2} \hat{\partial}_0 g_{jj} = \mathcal{O}(c^{-5}), \qquad (104)$$

<footnote n="26">²⁶ On top of the strong isotropy condition imposed on the metric, which fixes the space coordinate (Damour et al. 1991; Soffel 2003)—see Sect. 2.1.7.</footnote>
$$ \partial_j g_{ij} - \frac{1}{2} \partial_i (g_{jj} - g_{00}) = \mathcal{O}(c^{-4}). \qquad (105) $$

It implies that the field equations are modified, such that, instead of Eqs. (81)–(82), one has (Damour et al. 1991)

$$ \Delta w^{\text{SPN}} = -4\pi G\sigma + \mathcal{O}(c^{-4}), \qquad (106) $$

$$ \Delta w_i^{\text{SPN}} - \frac{1}{4} \partial_t \partial_i w^{\text{SPN}} = -4\pi G \sigma^i + \mathcal{O}(c^{-2}). \qquad (107) $$

The difference between the potentials in the harmonic and the SPN gauges reads (Soffel 2003)

$$ w^{\text{SPN}} = w + c^{-2} \partial_{tt}^2 \chi / 2 + \mathcal{O}(c^{-3}), \quad w_i^{\text{SPN}} = w_i - \partial_{tt}^2 \chi / 8, \qquad (108) $$

with

$$ \chi := G \int d^3 x' \sigma(t, x') \|x - x'\|. \qquad (109) $$

In terms of coordinates, the transformation between the two coordinate systems reads (Klerion and Soffel 2000; Will 2018b)

$$ t^{\text{SPN}} = t - \frac{1}{2} \partial_t \chi, \qquad (110) $$

$$ \mathbf{x}^{\text{SPN}} = \mathbf{x}. \qquad (111) $$

While both gauges can be used in principle, the harmonic gauge has been found to be much more convenient for deriving the reference frame theory in the solar system (Damour et al. 1991), such that harmonic coordinates have been used in the reference frame theory upon which the IAU recommendations have been built. A thorough comparison between the two coordinate systems in the full PPN framework can be found in Will (2018b).

## Full PPN metric in harmonic gauge

The full PPN framework considers all the possible different potentials that could contribute to the metric, assuming that the Newtonian potential is the leading potential—which may not be the case for a massive gravity theory for instance, see Sect. 4.4. Furthermore, the PPN framework can only deal with theories for which the Weak Equivalence Principle has been enforced by hand by requiring that the additional gravitational fields do not appear in the material part of the action. In the PPN framework, test particles indeed follow the geodesics of space-time, meaning that matter reacts to the space-time metric only. While this is necessary in order to have a theory for which all test particles fall alike by construction, one could question
$$ \partial_j g_{ij} - \frac{1}{2} \partial_i (g_{jj} - g_{00}) = \mathcal{O}(c^{-4}). \qquad (105) $$

It implies that the field equations are modified, such that, instead of Eqs. (81)–(82), one has (Damour et al. 1991)

$$ \Delta w^{\text{SPN}} = -4\pi G\sigma + \mathcal{O}(c^{-4}), \qquad (106) $$

$$ \Delta w_i^{\text{SPN}} - \frac{1}{4} \partial_t \partial_i w^{\text{SPN}} = -4\pi G \sigma^i + \mathcal{O}(c^{-2}). \qquad (107) $$

The difference between the potentials in the harmonic and the SPN gauges reads (Soffel 2003)

$$ w^{\text{SPN}} = w + c^{-2} \partial_{tt}^2 \chi / 2 + \mathcal{O}(c^{-3}), \quad w_i^{\text{SPN}} = w_i - \partial_{tt}^2 \chi / 8, \qquad (108) $$

with

$$ \chi := G \int d^3 x' \sigma(t, x') \|x - x'\|. \qquad (109) $$

In terms of coordinates, the transformation between the two coordinate systems reads (Klerion and Soffel 2000; Will 2018b)

$$ t^{\text{SPN}} = t - \frac{1}{2} \partial_t \chi, \qquad (110) $$

$$ \mathbf{x}^{\text{SPN}} = \mathbf{x}. \qquad (111) $$

While both gauges can be used in principle, the harmonic gauge has been found to be much more convenient for deriving the reference frame theory in the solar system (Damour et al. 1991), such that harmonic coordinates have been used in the reference frame theory upon which the IAU recommendations have been built. A thorough comparison between the two coordinate systems in the full PPN framework can be found in Will (2018b).

## Full PPN metric in harmonic gauge

The full PPN framework considers all the possible different potentials that could contribute to the metric, assuming that the Newtonian potential is the leading potential—which may not be the case for a massive gravity theory for instance, see Sect. 4.4. Furthermore, the PPN framework can only deal with theories for which the Weak Equivalence Principle has been enforced by hand by requiring that the additional gravitational fields do not appear in the material part of the action. In the PPN framework, test particles indeed follow the geodesics of space-time, meaning that matter reacts to the space-time metric only. While this is necessary in order to have a theory for which all test particles fall alike by construction, one could question
the theoretical relevance of such a restriction (Damour 2012)—see Weak Equivalence Principle in Sect. 4.2.

The full PPN metric in harmonic coordinates<fnref n="27" /> reads (Will 2018b):

$$g_{00} = -1 + 2U + [2\psi_{\text{Harm}} - 2\beta U^2] + \left(1 - \frac{1}{2}\alpha_1 + \alpha_2 + 2\xi\right)\ddot{X} + \Phi_{\text{Harm}}^{\text{PF}}, \quad (112)$$

$$g_{0j} = - \left[ 2(1 + \gamma) + \frac{1}{2}\alpha_1 \right] V_j - \frac{1}{4}\alpha_1 X_{,0j} + \Phi_{\text{Harm}}^{\text{PF}}, \quad (113)$$

$$g_{jk} = (1 + 2\gamma U)\delta_{jk} \quad (114)$$

with

$$\psi_{\text{Harm}} = \frac{1}{2}(2\gamma + 1 - 2\xi)\Phi_1 - (2\beta - 1 - \xi)\Phi_2 + \Phi_3 + (3\gamma - 2\xi)\Phi_4 \\ + \xi\Phi_6 - \xi\Phi_W, \quad (115)$$

$$\Phi_{\text{Harm}}^{\text{PF}} = -\alpha_1 w^2 U - \alpha_1 w^j V_j + \alpha_2 w^j w^k X_{jk} - \left(2\alpha_2 - \frac{1}{2}\alpha_1\right) w^j X_{,0j}, \quad (116)$$

$$\Phi_{\text{Harm}}^{\text{PF}} = -\frac{1}{2}\alpha_1 w^j U + \frac{1}{4}\alpha_1 w^k X_{jk}, \quad (117)$$

where $w$ is the coordinate velocity of the PPN coordinate system relative to the mean rest frame of the universe, and with the following metric potentials

$$U = \int \frac{\rho^{*}}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', V_j = \int \frac{\rho^{*}v_j}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', X = \int \rho^{*} \|\mathbf{x} - \mathbf{x}'\| d^3\mathbf{x}', \quad (118)$$

$$\Phi_1 = \int \frac{\rho^{*} v^2}{\|\mathbf{x} - \mathbf{x}'\|} d^3 \mathbf{x}', \quad \Phi_6 = \int \frac{\rho^{*} [\mathbf{v}' \cdot (\mathbf{x} - \mathbf{x}')]^2}{\|\mathbf{x} - \mathbf{x}'\|^3} d^3 \mathbf{x}', \quad (119)$$

$$\Phi_2 = \int \frac{\rho^{*} U'}{\|\mathbf{x} - \mathbf{x}'\|} d^3 \mathbf{x}', \quad \Phi_3 = \int \frac{\rho^{*} \Pi'}{\|\mathbf{x} - \mathbf{x}'\|} d^3 \mathbf{x}', \quad (120)$$

$$\Phi_4 = \int \frac{p'}{\|\mathbf{x} - \mathbf{x}'\|} d^3 \mathbf{x}', \quad (121)$$

$$\Phi_W = \int \rho^{*'} \rho^{*''} \frac{(\mathbf{x} - \mathbf{x}')}{\|\mathbf{x} - \mathbf{x}'\|^3} \cdot \left[ \frac{(\mathbf{x}' - \mathbf{x}'')}{\|\mathbf{x} - \mathbf{x}'\|} - \frac{(\mathbf{x} - \mathbf{x}'')}{\|\mathbf{x}' - \mathbf{x}''\|} \right] d^3 \mathbf{x}' d^3 \mathbf{x}'' \quad (122)$$

<footnote n="27">²⁷ Defined here as a “Nutku gauge”, that is $\partial_v[(1 + (1 + \gamma)U)\sqrt{-gg^{\mu\nu}}] = 0$, and not $\partial_v[\sqrt{-gg^{\mu\nu}}] = 0$— see discussion on the ambiguity of the definition of harmonic coordinates in Sect. 4.1.</footnote>
the theoretical relevance of such a restriction (Damour 2012)—see Weak Equivalence Principle in Sect. 4.2.

The full PPN metric in harmonic coordinates<fnref n="27" /> reads (Will 2018b):

$$g_{00} = -1 + 2U + [2\psi_{\text{Harm}} - 2\beta U^2] + \left(1 - \frac{1}{2}\alpha_1 + \alpha_2 + 2\xi\right)\ddot{X} + \Phi_{\text{Harm}}^{\text{PF}}, \quad (112)$$

$$g_{0j} = - \left[ 2(1 + \gamma) + \frac{1}{2}\alpha_1 \right] V_j - \frac{1}{4}\alpha_1 X_{,0j} + \Phi_{\text{Harm}}^{\text{PF}}, \quad (113)$$

$$g_{jk} = (1 + 2\gamma U)\delta_{jk} \quad (114)$$

with

$$\psi_{\text{Harm}} = \frac{1}{2}(2\gamma + 1 - 2\xi)\Phi_1 - (2\beta - 1 - \xi)\Phi_2 + \Phi_3 + (3\gamma - 2\xi)\Phi_4 + \xi\Phi_6 - \xi\Phi_W, \quad (115)$$

$$\Phi_{\text{Harm}}^{\text{PF}} = -\alpha_1 w^2 U - \alpha_1 w^j V_j + \alpha_2 w^j w^k X_{jk} - \left(2\alpha_2 - \frac{1}{2}\alpha_1\right) w^j X_{,0j}, \quad (116)$$

$$\Phi_{\text{Harm}}^{\text{PF}} = -\frac{1}{2}\alpha_1 w^j U + \frac{1}{4}\alpha_1 w^k X_{jk}, \quad (117)$$

where $w$ is the coordinate velocity of the PPN coordinate system relative to the mean rest frame of the universe, and with the following metric potentials

$$U = \int \frac{\rho^{*'}}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', V_j = \int \frac{\rho^{*'} v_j}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', X = \int \rho^{*'} \|\mathbf{x} - \mathbf{x}'\| d^3\mathbf{x}', \quad (118)$$

$$\Phi_1 = \int \frac{\rho^{*'} v^2}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', \quad \Phi_6 = \int \frac{\rho^{*'} [\mathbf{v}' \cdot (\mathbf{x} - \mathbf{x}')]^2}{\|\mathbf{x} - \mathbf{x}'\|^3} d^3\mathbf{x}', \quad (119)$$

$$\Phi_2 = \int \frac{\rho^{*'} U'}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', \quad \Phi_3 = \int \frac{\rho^{*'} \Pi'}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', \quad (120)$$

$$\Phi_4 = \int \frac{p'}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', \quad (121)$$

$$\Phi_W = \int \rho^{*'} \rho^{*''} \frac{(\mathbf{x} - \mathbf{x}')}{\|\mathbf{x} - \mathbf{x}'\|^3} \cdot \left[ \frac{(\mathbf{x}' - \mathbf{x}'')}{\|\mathbf{x} - \mathbf{x}'\|} - \frac{(\mathbf{x} - \mathbf{x}'')}{\|\mathbf{x}' - \mathbf{x}''\|} \right] d^3\mathbf{x}' d^3\mathbf{x}'' \quad (122)$$

<footnote n="27">²⁷ Defined here as a “Nutku gauge”, that is $\partial_v[(1 + (1 + \gamma)U)\sqrt{-gg^{\mu\nu}}] = 0$, and not $\partial_v[\sqrt{-gg^{\mu\nu}}] = 0$— see discussion on the ambiguity of the definition of harmonic coordinates in Sect. 4.1.</footnote>
where $\rho^* := \rho u^0 \sqrt{-g}$ is the “conserved” mass density$^{28}$—in the sense that it satisfies the Newtonian conservation equation $\partial_t \rho^* + \partial_j (\rho^* v^j) = 0$..$^{29}$ The whole set of PPN parameters is $\gamma, \beta, \xi, \alpha_1, \alpha_2, \alpha_3, \zeta_1, \zeta_2, \zeta_3, \zeta_4$. The parameter $\xi$ is non-zero in any theory of gravity that predicts preferred-location effects such as a galaxy-induced anisotropy in the local gravitational constant GL (also called “Whitehead” effects); $\alpha_1, \alpha_2, \alpha_3$ measure whether or not the theory predicts post-Newtonian preferred-frame effects; $\alpha_3, \zeta_1, \zeta_2, \zeta_3, \zeta_4$ measure whether or not the theory predicts violations of global conservation laws for total momentum (Will 2014).

Untill now, the full PPN framework has not been investigated with planetary ephemerides. First, one would have to derive the equation of motion in harmonic coordinates for point-mass interations that follows from the more general hydrodynamical situation depicted in Eqs. (112)–(122). Also, let us stress that the issue of the definition of the SSB has to be carefully checked, notably for the parameters that correspond to a violation of the conservation laws of momentum. Finally, the more the parameters, the more computationally demanding the study is, and one therefore has to think about the best way(s) to explore and constrain the landscape of parameters. Hence, several studies still seem to be needed before being able to constrain the full PPN framework with planetary ephemerides.

Nevertheless, measurements of PPN parameters other than $\beta$ and $\gamma$ such as $\alpha_1, \alpha_2, \zeta_1$, in the weak field regime, have been proposed early, by Anderson et al. (1978), and in the context of the mission BepiColombo by Milani et al. (2002) and De Marchi and Cascioli (2020), using Mercury orbit. As explained in Will (2014), the contribution of the preferred-frame parameters are included in the ratio between the planet mass and the Sun mass, leading to contributions to a maximum of $1 \times 10^{-4}$ for Jupiter and $1 \times 10^{-7}$ for Mercury. De Marchi and Cascioli (2020) expect a constraint for $\alpha_1$ at the level of $10^{-7}$ and up to $7 \times 10^{-8}$ for $\alpha_2$, with BepiColombo and future Venus missions. These estimations are obtained with covariance analysis of tracking observations of s/c orbiting planets. They are interesting as they estimate the capabilities of future missions for disentangling of the PPN contributions on the planetary motions.

In the strong field regime, constraints have been obtained for the parameters $\alpha_1, \alpha_3$ and $\zeta_2$ using pulsar timing (Stairs 2003).

## 4.2 Equivalence principle

While a deviation for particles with different compositions falling equally in a given gravitational field would indicate a violation of the Weak Equivalence Principle (WEP), a deviation for extended bodies with different gravitational self-energies would be a sign of a violation of the Gravitational Weak Equivalence Principle (GWEP), which is part of the Strong Equivalence Principle (SEP) as defined in Will (2018b). In all cases, with planetary ephemerides, one seeks to check the Universality of Free Fall (UFF)—that is, to check whether or not all bodies fall alike in a given

<footnote>$^{28}$ $\rho$ is the density of rest mass of a fluid element as measured in a local, freely falling, momentarily comoving frame. One has $\rho = \sigma + \mathcal{O}(c^{-2})$, where $\sigma$ is defined in Eq. (83).</footnote>

<footnote>$^{29}$ Note that one therefore has $X = \chi + \mathcal{O}(c^{-2})$.</footnote>
where $\rho^* := \rho u^0 \sqrt{-g}$ is the “conserved” mass density$^{28}$—in the sense that it satisfies the Newtonian conservation equation $\partial_t \rho^* + \partial_j (\rho^* v^j) = 0$..$^{29}$ The whole set of PPN parameters is $\gamma, \beta, \xi, \alpha_1, \alpha_2, \alpha_3, \zeta_1, \zeta_2, \zeta_3, \zeta_4$. The parameter $\xi$ is non-zero in any theory of gravity that predicts preferred-location effects such as a galaxy-induced anisotropy in the local gravitational constant GL (also called “Whitehead” effects); $\alpha_1, \alpha_2, \alpha_3$ measure whether or not the theory predicts post-Newtonian preferred-frame effects; $\alpha_3, \zeta_1, \zeta_2, \zeta_3, \zeta_4$ measure whether or not the theory predicts violations of global conservation laws for total momentum (Will 2014).

Untill now, the full PPN framework has not been investigated with planetary ephemerides. First, one would have to derive the equation of motion in harmonic coordinates for point-mass interations that follows from the more general hydrodynamical situation depicted in Eqs. (112)–(122). Also, let us stress that the issue of the definition of the SSB has to be carefully checked, notably for the parameters that correspond to a violation of the conservation laws of momentum. Finally, the more the parameters, the more computationally demanding the study is, and one therefore has to think about the best way(s) to explore and constrain the landscape of parameters. Hence, several studies still seem to be needed before being able to constrain the full PPN framework with planetary ephemerides.

Nevertheless, measurements of PPN parameters other than $\beta$ and $\gamma$ such as $\alpha_1, \alpha_2, \zeta_1$, in the weak field regime, have been proposed early, by Anderson et al. (1978), and in the context of the mission BepiColombo by Milani et al. (2002) and De Marchi and Cascioli (2020), using Mercury orbit. As explained in Will (2014), the contribution of the preferred-frame parameters are included in the ratio between the planet mass and the Sun mass, leading to contributions to a maximum of $1 \times 10^{-4}$ for Jupiter and $1 \times 10^{-7}$ for Mercury. De Marchi and Cascioli (2020) expect a constraint for $\alpha_1$ at the level of $10^{-7}$ and up to $7 \times 10^{-8}$ for $\alpha_2$, with BepiColombo and future Venus missions. These estimations are obtained with covariance analysis of tracking observations of s/c orbiting planets. They are interesting as they estimate the capabilities of future missions for disentangling of the PPN contributions on the planetary motions.

In the strong field regime, constraints have been obtained for the parameters $\alpha_1, \alpha_3$ and $\zeta_2$ using pulsar timing (Stairs 2003).

## 4.2 Equivalence principle

While a deviation for particles with different compositions falling equally in a given gravitational field would indicate a violation of the Weak Equivalence Principle (WEP), a deviation for extended bodies with different gravitational self-energies would be a sign of a violation of the Gravitational Weak Equivalence Principle (GWEP), which is part of the Strong Equivalence Principle (SEP) as defined in Will (2018b). In all cases, with planetary ephemerides, one seeks to check the Universality of Free Fall (UFF)—that is, to check whether or not all bodies fall alike in a given

<footnote>$^{28}$ $\rho$ is the density of rest mass of a fluid element as measured in a local, freely falling, momentarily comoving frame. One has $\rho = \sigma + \mathcal{O}(c^{-2})$, where $\sigma$ is defined in Eq. (83).</footnote>

<footnote>$^{29}$ Note that one therefore has $X = \chi + \mathcal{O}(c^{-2})$.</footnote>
gravitational potential. Such violations naturally occur in various modifications of general relativity, and in particular in theories with more than four dimensions like in string theories (Damour and Polyakov 1994; Damour et al. 2002), for which the constants of the standard model of particles turn out to be dynamical entities (Damour and Donoghue 2010; Damour 2012).

As we will see, the violations of both the WEP and GWEP have in common the fact that the equation of motion at the Newtonian level reads as follows

$$ \boldsymbol{a}_{\mathrm{T}}=-\sum_{A \neq T} \mu_{A}^{G} \boldsymbol{r}_{A T} \quad(123) $$

where $\delta_T$ and $\delta_{AT}$ are coefficients that depend on the composition of the bodies T and A, for the case of the violation of the WEP only. For the case of the violation of the GWEP only, $\delta_T$ depends on the ratio between the gravitational and inertial masses and $\delta_{AT} = 0$—such that

$$ \left(\frac{m^{G}}{m^{\prime}}\right)_{T}=1+\delta_{T}, \quad(124) $$

where $m^G$ and $m^\prime$ are the gravitational and inertial masses respectively. Note that $\mu_A^G$ is the gravitational parameter constructed with the gravitational mass $m_A^G$ such as $\mu_A^G = (1 + \delta_A)\mu_A^\prime$, where $\mu_A^\prime := Gm_A^\prime$ is the gravitational parameter obtained with the inertial mass $m_A^\prime$. Equation (123) is the most general equation that still satisfies Newton's third law of motion (Viswanathan et al. 2018). As we will see in Sect. 4.6, one can have both a violation of the WEP and the GWEP at the same time—that is, $\delta_T$ can depend both on the composition of the body T and on its internal gravitational energy.

The equivalence principle is currently tested at different scales: on the one hand on the laboratory scale for WEP with the torsion balance and MICROSCOPE, and on the other hand, the astronomical scale for GWEP with LLR<fnref n="30" /> and compact objects — at the level of about $10^{-13}$ with torsion balances (Adelberger et al. 2003a), $10^{-15}$ with the MICROSCOPE experiment (Touboul et al. 2022) and $10^{-14}$ with LLR (Williams et al. 2012; Viswanathan et al. 2018; Biskupek et al. 2021).

The latest LLR results (Biskupek et al. 2021) are more accurate that what can be achieved with the sensitivity of planetary ephemerides as discussed in Viswanathan et al. (2018). We give however a brief overview of the equations at stake below, mainly in the perspective of the BepiColombo mission that should lead to a significant improvement of the sensitivity of planetary ephemerides on the tests of the equivalence principle.

Finally, it is interesting to note that constraints on the violation of SEP have been obtained in strong field regime, mainly using the chronometry of millisecond pulsars. For a deeper dive into this captivating subject, we direct the reader to Kramer et al. (2021), Freire et al. (2012), and Voisin et al. (2020).

<footnote n="30">³⁰ LLR tests also include a WEP contribution as the Earth and Moon have different compositions (Adelberger 2001).</footnote>
gravitational potential. Such violations naturally occur in various modifications of general relativity, and in particular in theories with more than four dimensions like in string theories (Damour and Polyakov 1994; Damour et al. 2002), for which the constants of the standard model of particles turn out to be dynamical entities (Damour and Donoghue 2010; Damour 2012).

As we will see, the violations of both the WEP and GWEP have in common the fact that the equation of motion at the Newtonian level reads as follows

$$ \mathbf{a}_T = -\sum_{A \neq T} \mu_A^G \frac{\mathbf{r}_{AT}}{r_{AT}^3} (1 + \delta_T + \delta_{AT}), \qquad (123) $$

where $\delta_T$ and $\delta_{AT}$ are coefficients that depend on the composition of the bodies T and A, for the case of the violation of the WEP only. For the case of the violation of the GWEP only, $\delta_T$ depends on the ratio between the gravitational and inertial masses and $\delta_{AT} = 0$—such that

$$ \left( \frac{m^G}{m^I} \right)_T = 1 + \delta_T, \qquad (124) $$

where $m^G$ and $m^I$ are the gravitational and inertial masses respectively. Note that $\mu_A^G$ is the gravitational parameter constructed with the gravitational mass $m_A^G$ such as $\mu_A^G = (1 + \delta_A)\mu_A^I$, where $\mu_A^I := Gm_A^I$ is the gravitational parameter obtained with the inertial mass $m_A^I$. Equation (123) is the most general equation that still satisfies Newton's third law of motion (Viswanathan et al. 2018). As we will see in Sect. 4.6, one can have both a violation of the WEP and the GWEP at the same time—that is, $\delta_T$ can depend both on the composition of the body T and on its internal gravitational energy.

The equivalence principle is currently tested at different scales: on the one hand on the laboratory scale for WEP with the torsion balance and MICROSCOPE, and on the other hand, the astronomical scale for GWEP with LLR<fnref n="30" /> and compact objects — at the level of about $10^{-13}$ with torsion balances (Adelberger et al. 2003a), $10^{-15}$ with the MICROSCOPE experiment (Touboul et al. 2022) and $10^{-14}$ with LLR (Williams et al. 2012; Viswanathan et al. 2018; Biskupek et al. 2021).

The latest LLR results (Biskupek et al. 2021) are more accurate that what can be achieved with the sensitivity of planetary ephemerides as discussed in Viswanathan et al. (2018). We give however a brief overview of the equations at stake below, mainly in the perspective of the BepiColombo mission that should lead to a significant improvement of the sensitivity of planetary ephemerides on the tests of the equivalence principle.

Finally, it is interesting to note that constraints on the violation of SEP have been obtained in strong field regime, mainly using the chronometry of millisecond pulsars. For a deeper dive into this captivating subject, we direct the reader to Kramer et al. (2021), Freire et al. (2012), and Voisin et al. (2020).

<footnote n="30">³⁰ LLR tests also include a WEP contribution as the Earth and Moon have different compositions (Adelberger 2001).</footnote>
### 4.2.1 The case of metric theories

While metric theories satisfy the WEP by design (see Sect. 4.1), they do not satisfy the SEP in general—in the sense that bodies with different gravitational self-energies would not follow the same trajectories in general (Nordtvedt 1968a, b). This notably induces that for extended bodies with enough gravitational self-energy such as planets, the ratio between their gravitational and the inertial masses is no longer one as in general relativity, but instead depends on the body self-gravitational energy that can be approximated for uniformly distributed spherical bodies as follows

$$ \left(\frac{m^G}{m^I}\right)_T^{SEP} = 1 + \eta \times \frac{3 G m_T^G}{5 c^2 R_T}, \qquad (125) $$

where $R_T$ is the radius of the astronomical body T and $m_T^G$ is its gravitational mass. The parameter $\eta$, known at the Nordtvedt parameter, quantifies the possible violation of the SEP and can be fitted together with the rest of the planetary and lunar ephemeris parameters. In the most generic case, $\eta$ can be seen as an independent parameter and be introduced in Eq. (125), common to all the bodies. However, in the case of the metric theories, it can be related to the PPN parameters such as (Will 2018b):

$$ \eta = 4(\beta - 1) - (\gamma - 1) - \alpha_1 - \frac{2}{3}\alpha_2. \qquad (126) $$

Although note that only the cases for which $\alpha_1 = 0$ and $\alpha_2 = 0$ have been considered in planetary ephemerides so far—see Sect. 4.1.4 for a discussion on the extended parametrized post-Newtonian framework. In principle, the introduction of Eq. (126) in the global adjustment of planetary ephemerides should increase the correlations between $\beta$, $\gamma$, $\eta$ and the other ephemeris parameters, including $J_2^\circ$ (Ashby et al. 2007). However, simulations for the MORE BepiColombo experiment indicate a significant improvement of the $\beta$, $J_2^\circ$ and $\eta$ correlation (from 0.9 to 0.3) thanks to the high accuracy of the Mercury-Earth range measurements that should be obtained (few tens of centimeters) by the mission (Milani et al. 2002; De Marchi et al. 2016; Imperi et al. 2018; De Marchi and Cascioli 2020). Finally, recent attempt using Messenger tracking data leads to a measurement of $\eta$ at the level of $7 \times 10^{-5}$ but this estimation has been done in considering only the Mercury orbit (Genova et al. 2018). Even in these very favorable conditions, it gives a slightly less restrictive limit than the latest LLR evaluation from Biskupek et al. (2021) with a deduced $\eta$ value of $5 \times 10^{-5}$.

### 4.2.2 The case of non-metric theories

As discussed in Sect. 4.1, metric theories satisfy the WEP because the additional gravitational fields (beyond the metric) do not couple directly to matter. However, one could question the reason why additional gravitational fields would not couple to matter. Indeed, at the fundamental level, there does not seem to be any reason to expect that outcome (Damour 2012). Moreover, extra-gravitational fields with non-
### 4.2.1 The case of metric theories

While metric theories satisfy the WEP by design (see Sect. 4.1), they do not satisfy the SEP in general—in the sense that bodies with different gravitational self-energies would not follow the same trajectories in general (Nordtvedt 1968a, b). This notably induces that for extended bodies with enough gravitational self-energy such as planets, the ratio between their gravitational and the inertial masses is no longer one as in general relativity, but instead depends on the body self-gravitational energy that can be approximated for uniformly distributed spherical bodies as follows

$$ \left(\frac{m^G}{m^I}\right)_T^{SEP} = 1 + \eta \times \frac{3 G m_T^G}{5 c^2 R_T}, \qquad (125) $$

where $R_T$ is the radius of the astronomical body T and $m_T^G$ is its gravitational mass. The parameter $\eta$, known at the Nordtvedt parameter, quantifies the possible violation of the SEP and can be fitted together with the rest of the planetary and lunar ephemeris parameters. In the most generic case, $\eta$ can be seen as an independent parameter and be introduced in Eq. (125), common to all the bodies. However, in the case of the metric theories, it can be related to the PPN parameters such as (Will 2018b):

$$ \eta = 4(\beta - 1) - (\gamma - 1) - \alpha_1 - \frac{2}{3}\alpha_2. \qquad (126) $$

Although note that only the cases for which $\alpha_1 = 0$ and $\alpha_2 = 0$ have been considered in planetary ephemerides so far—see Sect. 4.1.4 for a discussion on the extended parametrized post-Newtonian framework. In principle, the introduction of Eq. (126) in the global adjustment of planetary ephemerides should increase the correlations between $\beta$, $\gamma$, $\eta$ and the other ephemeris parameters, including $J_2^\circ$ (Ashby et al. 2007). However, simulations for the MORE BepiColombo experiment indicate a significant improvement of the $\beta$, $J_2^\circ$ and $\eta$ correlation (from 0.9 to 0.3) thanks to the high accuracy of the Mercury-Earth range measurements that should be obtained (few tens of centimeters) by the mission (Milani et al. 2002; De Marchi et al. 2016; Imperi et al. 2018; De Marchi and Cascioli 2020). Finally, recent attempt using Messenger tracking data leads to a measurement of $\eta$ at the level of $7 \times 10^{-5}$ but this estimation has been done in considering only the Mercury orbit (Genova et al. 2018). Even in these very favorable conditions, it gives a slightly less restrictive limit than the latest LLR evaluation from Biskupek et al. (2021) with a deduced $\eta$ value of $5 \times 10^{-5}$.

### 4.2.2 The case of non-metric theories

As discussed in Sect. 4.1, metric theories satisfy the WEP because the additional gravitational fields (beyond the metric) do not couple directly to matter. However, one could question the reason why additional gravitational fields would not couple to matter. Indeed, at the fundamental level, there does not seem to be any reason to expect that outcome (Damour 2012). Moreover, extra-gravitational fields with non-
minimal couplings to matter are ubiquitous in various attempts to unify the whole set of fundamental interactions in physics—such as in Kaluza-Klein or superstring theories (Damour 2012), for instance. Finally, such non-minimal couplings can also appear from radiative quantum corrections of theories that satisfy the WEP at the classical level (Armendariz-Picon and Penco 2012)—whereas the WEP in general relativity is immune to quantum corrections due to the specific symmetries of general relativity (Armendariz-Picon and Penco 2012). Therefore, it seems that beyond general relativity, violations of the WEP should be expected at some level—either at the classical level already, or from quantum corrections to the classical action.

### Decoupling mechanisms

Since the equivalence principle is observed to be satisfied with an ever increasing level of accuracy, non-metric theories need to possess mechanisms that are able to hide the effects of the additional gravitational fields in the solar system. They are often called “decoupling” or “screening” mechanisms.<fnref n="31" /> They can take their roots from the fact that the additional fields are simply massive enough to propagate over very short distances only—such as what happens with flux compactification in string theory (Douglas and Kachru 2007)—or from the suppression of the coupling during the evolution of the universe—such as the Damour–Nordtvedt decoupling mechanism (Damour and Nordtvedt 1993; Damour and Polyakov 1994)—or from the dependence of the effective mass, or coupling, of the additional gravitational field to the local density—such as the chameleon and symmetron mechanisms (Khoury and Weltman 2004; Hinterbichler and Khoury 2010)—or from non-perturbative effects— such as the Vainstein mechanism (Vainshtein 1972; Babichev and Deffayet 2013)— or simply from exact cancellations in the field equations—such as the intrinsic decoupling (Minazzoli and Hees 2013, 2016) that notably appears in entangled relativity (Arruga et al. 2021).

### 4.2.3 The equivalence principle in a planet-satellite system

Considering the generic Eq. (123), we can deduce for the couple of the bodies T and A that the difference between their accelerations toward the Sun S reads

$$ \Delta \mathbf{a} \equiv \mathbf{a}_A - \mathbf{a}_T = -\frac{G \bar{m}_{AT}}{r_{TA}^3} \mathbf{r}_{TA} + G m_S^G \left[ \frac{\mathbf{r}_{ST}}{r_{ST}^3} - \frac{\mathbf{r}_{SA}}{r_{SA}^3} \right] + G m_S^G \left[ \frac{\mathbf{r}_{ST}}{r_{ST}^3} (\delta_T + \delta_{ST}) - \frac{\mathbf{r}_{SA}}{r_{SA}^3} (\delta_A + \delta_{SA}) \right], \quad (127) $$

with $\bar{m}_{AT} \equiv m_A^G + m_T^G + (\delta_T + \delta_{TA})m_A^G + (\delta_A + \delta_{TA})m_T^G$. For gravitationally bounded systems like the Earth–Moon system, one can further approximate both distances appearing in this last term as being approximately equal, such that the relevant term in order to check a potential violation of the EP through a violation of the universality of free fall (UFF) is

<footnote n="31">³¹ Note that decoupling and screening mechanisms have been proposed first and foremost for theories that do satisfy the WEP, often in order to be able to explain the acceleration of the expansion of the universe without impacting the well-constrained solar system phenomenology.</footnote>
minimal couplings to matter are ubiquitous in various attempts to unify the whole set of fundamental interactions in physics—such as in Kaluza-Klein or superstring theories (Damour 2012), for instance. Finally, such non-minimal couplings can also appear from radiative quantum corrections of theories that satisfy the WEP at the classical level (Armendariz-Picon and Penco 2012)—whereas the WEP in general relativity is immune to quantum corrections due to the specific symmetries of general relativity (Armendariz-Picon and Penco 2012). Therefore, it seems that beyond general relativity, violations of the WEP should be expected at some level—either at the classical level already, or from quantum corrections to the classical action.

### Decoupling mechanisms

Since the equivalence principle is observed to be satisfied with an ever increasing level of accuracy, non-metric theories need to possess mechanisms that are able to hide the effects of the additional gravitational fields in the solar system. They are often called “decoupling” or “screening” mechanisms.<fnref n="31" /> They can take their roots from the fact that the additional fields are simply massive enough to propagate over very short distances only—such as what happens with flux compactification in string theory (Douglas and Kachru 2007)—or from the suppression of the coupling during the evolution of the universe—such as the Damour–Nordtvedt decoupling mechanism (Damour and Nordtvedt 1993; Damour and Polyakov 1994)—or from the dependence of the effective mass, or coupling, of the additional gravitational field to the local density—such as the chameleon and symmetron mechanisms (Khoury and Weltman 2004; Hinterbichler and Khoury 2010)—or from non-perturbative effects— such as the Vainstein mechanism (Vainshtein 1972; Babichev and Deffayet 2013)— or simply from exact cancellations in the field equations—such as the intrinsic decoupling (Minazzoli and Hees 2013, 2016) that notably appears in entangled relativity (Arruga et al. 2021).

### 4.2.3 The equivalence principle in a planet-satellite system

Considering the generic Eq. (123), we can deduce for the couple of the bodies T and A that the difference between their accelerations toward the Sun S reads

$$ \Delta \mathbf{a} \equiv \mathbf{a}_A - \mathbf{a}_T = - \frac{G \bar{m}_{AT}}{r_{TA}^3} \mathbf{r}_{TA} + G m_S^G \left[ \frac{\mathbf{r}_{ST}}{r_{ST}^3} - \frac{\mathbf{r}_{SA}}{r_{SA}^3} \right] + G m_S^G \left[ \frac{\mathbf{r}_{ST}}{r_{ST}^3} (\delta_T + \delta_{ST}) - \frac{\mathbf{r}_{SA}}{r_{SA}^3} (\delta_A + \delta_{SA}) \right], \quad (127) $$

with $\bar{m}_{AT} \equiv m_A^G + m_T^G + (\delta_T + \delta_{TA})m_A^G + (\delta_A + \delta_{TA})m_T^G$. For gravitationally bounded systems like the Earth–Moon system, one can further approximate both distances appearing in this last term as being approximately equal, such that the relevant term in order to check a potential violation of the EP through a violation of the universality of free fall (UFF) is

<footnote n="31">³¹ Note that decoupling and screening mechanisms have been proposed first and foremost for theories that do satisfy the WEP, often in order to be able to explain the acceleration of the expansion of the universe without impacting the well-constrained solar system phenomenology.</footnote>
$$ \Delta \mathbf{a}^{\overline{\text{UFF}}} \approx \mathbf{a}_E \Delta \text{ESM} \qquad (128) $$

where E, S and M stand for the Earth, the Sun and the Moon respectively, and with

$$ \Delta_{\text{ESM}} = [(\delta_E + \delta_{\text{SE}}) - (\delta_M + \delta_{\text{SM}})]. \qquad (129) $$

Equation (128) was introduced in Williams et al. (2012) for theories that violate the GWEP only, for which $\Delta_{\text{ESM}} = \Delta_{\text{EM}}$ such that

$$ \Delta_{\text{EM}} = \left[ \left( \frac{m^G}{m'} \right)_E - \left( \frac{m^G}{m'} \right)_M \right]. \qquad (130) $$

In that case, Eq. (128) becomes

$$ \Delta \mathbf{a}^{\overline{\text{UFF}}} \approx \mathbf{a}_E \Delta_{\text{EM}}, \qquad (131) $$

### 4.2.4 The equivalence principle at planetary scale

At planetary scale, Eqs. (128) or (131) are not valid anymore. Indeed, planets are at different distances from the Sun and EP tests cannot be limited to the differences between two accelerations. In that case, one has to use the more general Eq. (123) instead.

A confusion can occur when defining the SSB because the gravitational parameters that appear in the equations of motion are constructed with the gravitational masses; whereas the gravitational parameters that are used to define the SSB should be constructed with the inertial masses, as pointed out in Milani et al. (2002) and Genova et al. (2018)—see also Sect. 4.6. It was an issue in particular for previous versions of planetary ephemerides, such as Folkner et al. (2008), where the SSB positions and velocities were estimated at each step of integration with the gravitational parameters appearing in the equations of motion.

However, as explained in Sect. 3.2.4, in planetary ephemerides such as INPOP08 (Fienga et al. 2008) and the followings, the SSB position and velocity are estimated once, before integrating the motion of planets, using planetary initial conditions, inertial masses and Eqs. (89) and (88). As explained in Sect. 4.6, Bernus et al. (2022) have checked with a specific theory that violates both the WEP and the GWEP that those equations indeed remain the ones that characterize the SSB at the required level of accuracy for planetary ephemerides—provided that inertial masses are used for the definitions and not gravitational masses—as originally suggested in Milani et al. (2002) and Genova et al. (2018).

### 4.2.5 Comment on the mass in the Shapiro delay equation

As soon as the SEP is violated, there is a potential ambiguity about the gravitational parameter that appears in the Shapiro delay Eq. (97). Indeed, the gravitational
$$ \Delta \mathbf{a}^{\overline{\text{UFF}}} \approx \mathbf{a}_E \Delta \text{ESM} \qquad (128) $$

where E, S and M stand for the Earth, the Sun and the Moon respectively, and with

$$ \Delta_{\text{ESM}} = [(\delta_E + \delta_{\text{SE}}) - (\delta_M + \delta_{\text{SM}})]. \qquad (129) $$

Equation (128) was introduced in Williams et al. (2012) for theories that violate the GWEP only, for which $\Delta_{\text{ESM}} = \Delta_{\text{EM}}$ such that

$$ \Delta_{\text{EM}} = \left[ \left( \frac{m^G}{m'} \right)_E - \left( \frac{m^G}{m'} \right)_M \right]. \qquad (130) $$

In that case, Eq. (128) becomes

$$ \Delta \mathbf{a}^{\overline{\text{UFF}}} \approx \mathbf{a}_E \Delta_{\text{EM}}, \qquad (131) $$

### 4.2.4 The equivalence principle at planetary scale

At planetary scale, Eqs. (128) or (131) are not valid anymore. Indeed, planets are at different distances from the Sun and EP tests cannot be limited to the differences between two accelerations. In that case, one has to use the more general Eq. (123) instead.

A confusion can occur when defining the SSB because the gravitational parameters that appear in the equations of motion are constructed with the gravitational masses; whereas the gravitational parameters that are used to define the SSB should be constructed with the inertial masses, as pointed out in Milani et al. (2002) and Genova et al. (2018)—see also Sect. 4.6. It was an issue in particular for previous versions of planetary ephemerides, such as Folkner et al. (2008), where the SSB positions and velocities were estimated at each step of integration with the gravitational parameters appearing in the equations of motion.

However, as explained in Sect. 3.2.4, in planetary ephemerides such as INPOP08 (Fienga et al. 2008) and the followings, the SSB position and velocity are estimated once, before integrating the motion of planets, using planetary initial conditions, inertial masses and Eqs. (89) and (88). As explained in Sect. 4.6, Bernus et al. (2022) have checked with a specific theory that violates both the WEP and the GWEP that those equations indeed remain the ones that characterize the SSB at the required level of accuracy for planetary ephemerides—provided that inertial masses are used for the definitions and not gravitational masses—as originally suggested in Milani et al. (2002) and Genova et al. (2018).

### 4.2.5 Comment on the mass in the Shapiro delay equation

As soon as the SEP is violated, there is a potential ambiguity about the gravitational parameter that appears in the Shapiro delay Eq. (97). Indeed, the gravitational
parameter in the Shapiro delay could either be constructed on the inertial mass, or on the gravitational mass. As we will see in Sect. 4.6 Eqs. (157)–(158), one can derive from first principles that the gravitational parameter that appears in the Shapiro delay is the one built upon the inertial mass, that is $\mu_A^I = Gm_A^I = (1 - \delta_A)Gm_A^G$, with:

$$c(t_r - t_e) = R + \sum_A (\gamma + 1) \frac{\mu_A^I}{c^2} \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}, \quad (132)$$

$$= R + \sum_A (\gamma + 1) \frac{\mu_A^G}{c^2} \left( 1 - \eta \times \frac{3}{5} \frac{\mu_A^G}{c^2 R_A} \right) \times \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}. \quad (133)$$

There is a somewhat simple explanation to that. One must recall that the original Nordtvedt effect was derived for the sum of point particles of (inertial) mass $m$ that are gravitationally bounded together. The metric is bent by those masses $m$, and the sum of all the masses corresponds to the inertial mass of the overall body that is composed from all those masses, to lowest order (Nordtvedt 1968b). Since the Shapiro delay simply derives from the null condition $ds^2 = 0$ for light, the Shapiro delay must depend on this sum of inertial masses. What the Nordtvedt effect tells additionally, is that the center of mass of this collection of masses does not follow the same trajectory as would an inertial mass, but follows the trajectory given by the following equation instead

$$\mathbf{a}_T = - \sum_{A \neq T} \mu_A^G \frac{\mathbf{r}_{AT}}{r_{AT}^3} \left( 1 + \eta \times \frac{3}{5} \frac{\mu_T^G}{c^2 R_T} \right). \quad (134)$$

### 4.2.6 The issue of time in non-metric theories

The issue with theories that violate the WEP is that they usually also violate the local position invariance (Uzan 2011), such that the behaviour of clocks not only may be modified with respect to general relativity depending on the position in the gravitational field, but also may depend on the clock composition. Therefore, one may expect that this should complicate the tests of such theories with planetary ephemerides, given that positions of astronomical bodies are given in one of the coordinate times recommended by the IAU, and which are perfectly independent of the composition of the clocks that has been used. Fortunately, the local position invariance is tested to a very high level of accuracy on Earth with various atomic clocks (Rosenband et al. 2008; Guéna et al. 2012; Leefer et al. 2013; Godun et al. 2014; Huntemann et al. 2014). Indeed, thanks to the eccentricity of the Earth’s orbit, one can probe a potential variation of the ratio between the frequencies of atomic
parameter in the Shapiro delay could either be constructed on the inertial mass, or on the gravitational mass. As we will see in Sect. 4.6 Eqs. (157)–(158), one can derive from first principles that the gravitational parameter that appears in the Shapiro delay is the one built upon the inertial mass, that is $\mu_A^I = Gm_A^I = (1 - \delta_A)Gm_A^G$, with:

$$c(t_r - t_e) = R + \sum_A (\gamma + 1) \frac{\mu_A^I}{c^2} \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}, \quad (132)$$

$$= R + \sum_A (\gamma + 1) \frac{\mu_A^G}{c^2} \left( 1 - \eta \times \frac{3}{5} \frac{\mu_A^G}{c^2 R_A} \right) \times \ln \frac{\mathbf{n} \cdot \mathbf{r}_{rA} + r_{rA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}{\mathbf{n} \cdot \mathbf{r}_{eA} + r_{eA} + 2 \frac{(1+\gamma)\mu_A}{c^2}}. \quad (133)$$

There is a somewhat simple explanation to that. One must recall that the original Nordtvedt effect was derived for the sum of point particles of (inertial) mass $m$ that are gravitationally bounded together. The metric is bent by those masses $m$, and the sum of all the masses corresponds to the inertial mass of the overall body that is composed from all those masses, to lowest order (Nordtvedt 1968b). Since the Shapiro delay simply derives from the null condition $ds^2 = 0$ for light, the Shapiro delay must depend on this sum of inertial masses. What the Nordtvedt effect tells additionally, is that the center of mass of this collection of masses does not follow the same trajectory as would an inertial mass, but follows the trajectory given by the following equation instead

$$\mathbf{a}_T = - \sum_{A \neq T} \mu_A^G \frac{\mathbf{r}_{AT}}{r_{AT}^3} \left( 1 + \eta \times \frac{3}{5} \frac{\mu_T^G}{c^2 R_T} \right). \quad (134)$$

### 4.2.6 The issue of time in non-metric theories

The issue with theories that violate the WEP is that they usually also violate the local position invariance (Uzan 2011), such that the behaviour of clocks not only may be modified with respect to general relativity depending on the position in the gravitational field, but also may depend on the clock composition. Therefore, one may expect that this should complicate the tests of such theories with planetary ephemerides, given that positions of astronomical bodies are given in one of the coordinate times recommended by the IAU, and which are perfectly independent of the composition of the clocks that has been used. Fortunately, the local position invariance is tested to a very high level of accuracy on Earth with various atomic clocks (Rosenband et al. 2008; Guéna et al. 2012; Leefer et al. 2013; Godun et al. 2014; Huntemann et al. 2014). Indeed, thanks to the eccentricity of the Earth’s orbit, one can probe a potential variation of the ratio between the frequencies of atomic
clock with different compositions depending on the variation of the local value of the Sun's gravitational field. Because the local position invariance is tested with atomic clocks at a level of accuracy that is far beyond what is used in space navigation, one can neglect this potential effect in planetary ephemerides.

## 4.3 Variation of the gravitational constant G

With respect to the subject of the potential variation of constants (Uzan 2011), planetary and lunar ephemerides have been mostly used for testing an hypothetical variation of Newton's constant G, or a related value, the Astronomical Unit (Hellings et al. 1983; Anderson et al. 1987, 1989; Pitjeva 1993).

Before 2012 and the decision by the International Astronomical Union (Capitaine et al. 2012) to redefine the Astronomical Unit (AU) as a constant with a fixed value, the AU was indeed part of the parameters estimated during the construction of planetary ephemerides. The Gauss constant k was fixed (Standish 2001) and the AU was estimated from this equation

$$GM_{\odot} = (AU)^3 k^2 / D^2 \qquad (135)$$

where D is the length of the day and $GM_{\odot}$ is the gravitational parameter of the Sun. In this context, AU values and its hypothetical time variations have been estimated with planetary fit by Krasinsky and Brumberg (2004). In the same manner, an estimation of $\dot{G}$ was also obtained, together with the value of AU. A first test of such a variation has been initiated by Hellings et al. (1983), using the data obtained with the Viking mission on Mars. On a regular basis, the planetary ephemerides were then published with updates for $\dot{G}$ based on a global adjustment including also AU and PPN parameters (Pitjeva 1993; Standish 2001; Pitjeva 2005b). Results are gathered in Table 6.

After 2012, the AU has been fixed to 149597870700 m, closing the door to the estimation of dAU/dt, but leading to a more consistent determination of the gravitational parameter of the Sun, $GM_{\odot}$ noted $\mu_{\odot}$, now fitted in planetary ephemerides. As it is reviewed in Uzan (2011), some theories propose that the gravitational constant can vary with time. Thus, an hypothetical variation of G in time has been introduced in the planetary ephemerides, but it is not straightforward since it affects both the planetary and the Sun gravitational parameters $\mu$.

For the planets, because they do not lose or gain mass, the time variation of G induces a direct variation of the gravitational parameter $\mu_P$ such $\dot{\mu}_P / \mu_P = \dot{G} / G$. For the Sun, its intrinsic mass loss $\dot{M}_{\odot}$ has also to be considered. As a consequence, the following equations are added to the system of equations integrated numerically during the planetary ephemerides construction

$$\frac{\dot{\mu}_{\odot}}{\mu_{\odot}} = \frac{\dot{G}}{G} + \frac{\dot{M}_{\odot}}{M_{\odot}}, \qquad (136)$$

where $\mu_{\odot}$ is the gravitational parameter of the Sun and
clock with different compositions depending on the variation of the local value of the Sun's gravitational field. Because the local position invariance is tested with atomic clocks at a level of accuracy that is far beyond what is used in space navigation, one can neglect this potential effect in planetary ephemerides.

## 4.3 Variation of the gravitational constant G

With respect to the subject of the potential variation of constants (Uzan 2011), planetary and lunar ephemerides have been mostly used for testing an hypothetical variation of Newton's constant G, or a related value, the Astronomical Unit (Hellings et al. 1983; Anderson et al. 1987, 1989; Pitjeva 1993).

Before 2012 and the decision by the International Astronomical Union (Capitaine et al. 2012) to redefine the Astronomical Unit (AU) as a constant with a fixed value, the AU was indeed part of the parameters estimated during the construction of planetary ephemerides. The Gauss constant k was fixed (Standish 2001) and the AU was estimated from this equation

$$GM_{\odot} = (AU)^3 k^2 / D^2 \qquad (135)$$

where D is the length of the day and $GM_{\odot}$ is the gravitational parameter of the Sun. In this context, AU values and its hypothetical time variations have been estimated with planetary fit by Krasinsky and Brumberg (2004). In the same manner, an estimation of $\dot{G}$ was also obtained, together with the value of AU. A first test of such a variation has been initiated by Hellings et al. (1983), using the data obtained with the Viking mission on Mars. On a regular basis, the planetary ephemerides were then published with updates for $\dot{G}$ based on a global adjustment including also AU and PPN parameters (Pitjeva 1993; Standish 2001; Pitjeva 2005b). Results are gathered in Table 6.

After 2012, the AU has been fixed to 149597870700 m, closing the door to the estimation of dAU/dt, but leading to a more consistent determination of the gravitational parameter of the Sun, $GM_{\odot}$ noted $\mu_{\odot}$, now fitted in planetary ephemerides. As it is reviewed in Uzan (2011), some theories propose that the gravitational constant can vary with time. Thus, an hypothetical variation of G in time has been introduced in the planetary ephemerides, but it is not straightforward since it affects both the planetary and the Sun gravitational parameters $\mu$.

For the planets, because they do not lose or gain mass, the time variation of G induces a direct variation of the gravitational parameter $\mu_P$ such $\dot{\mu}_P / \mu_P = \dot{G} / G$. For the Sun, its intrinsic mass loss $\dot{M}_{\odot}$ has also to be considered. As a consequence, the following equations are added to the system of equations integrated numerically during the planetary ephemerides construction

$$\frac{\dot{\mu}_{\odot}}{\mu_{\odot}} = \frac{\dot{G}}{G} + \frac{\dot{M}_{\odot}}{M_{\odot}}, \qquad (136)$$

where $\mu_{\odot}$ is the gravitational parameter of the Sun and
<table><caption>**Table 6** Values of $\dot{\mu}_{\odot}/\mu_{\odot}$ found in the literature deduced from planetary ephemerides. The Column 3 indicates the type of method used for the estimation: either the full fit (FF), the partial fit (PF) in or the random walk exploration (RW). See text for details. For Fienga et al. (2021a), PF<sup>**</sup> indicates that the oblateness of the Sun is strongly constrained by helio-seismological value. So this solution cannot be seen as a free fit but as a partial fit. Analysis of the MESSENGER mission is labeled MSG in Column 3 using only the Mercury orbit. The $\dot{G}/G$ values are deduced with $\dot{\mu}_{\odot}/\mu_{\odot}$ and $\dot{M}_{\odot}/M_{\odot}$ given in Column 5. The value of $\dot{G}/G$ deduced from Pitjeva and Pitjev (2013) is not the one published by the authors but was obtained in using the same value of $\dot{M}_{\odot}/M_{\odot}$ as for Konopliv et al. (2011)</caption><thead><tr><th>References</th><th>PE</th><th>Method</th><th>μ̇<sub>☉</sub>/μ<sub>☉</sub> × 10<sup>14</sup><br>yr<sup>-1</sup></th><th>Ṁ<sub>☉</sub>/M<sub>☉</sub> × 10<sup>14</sup><br>yr<sup>-1</sup></th><th>GW/G × 10<sup>14</sup><br>yr<sup>-1</sup></th></tr></thead><tbody><tr><td colspan="6"><em>Before AU fixed</em></td></tr><tr><td>Pitjeva (1993)</td><td>EPM1988</td><td>FF</td><td>NA</td><td>NA</td><td>470 ± 470</td></tr><tr><td>Pitjeva (2001)</td><td>EPM2000</td><td>FF</td><td>NA</td><td>NA</td><td>4 ± 8</td></tr><tr><td>Standish (2001)</td><td>DE405</td><td>FF</td><td>NA</td><td>NA</td><td>1 ± 8</td></tr><tr><td>Pitjeva (2005b)</td><td>EPM2004</td><td>FF</td><td>NA</td><td>NA</td><td>-2 ± 5</td></tr><tr><td colspan="6"><em>After AU fixed</em></td></tr><tr><td>Konopliv et al. (2011)</td><td>DE421</td><td>PF</td><td>1 ± 16</td><td>-9.2 ± 6.1</td><td>10.2 ± 22.1*</td></tr><tr><td>Pitjeva and Pitjev (2012)</td><td>EPM2010</td><td>FF</td><td>-5 ± 4</td><td>-6.7 ± 3.1</td><td>1.65 ± 8.77</td></tr><tr><td>Pitjeva and Pitjev (2013)</td><td>EPM2011</td><td>FF</td><td>-6.3 ± 6.4</td><td>-6.7 ± 3.1</td><td>0.4 ± 11.1</td></tr><tr><td rowspan="2">Fienga et al. (2015)</td><td rowspan="2">INPOP15a</td><td>FF</td><td>-5.0 ± 2.9</td><td>-9.2 ± 6.1</td><td>4.2 ± 9.0</td></tr><tr><td>RW</td><td>-4.3 ± 7.4</td><td>-9.2 ± 6.1</td><td>4.9 ± 13.5</td></tr><tr><td>Genova et al. (2018)</td><td>MSG</td><td></td><td>-6.130 ± 1.47</td><td>-10 ± 1</td><td>4 ± 7.5</td></tr><tr><td>Pitjeva et al. (2021)</td><td>EPM2021</td><td>FF</td><td>-10.2 ± 1.4</td><td>-11.05 ± 2.35</td><td>0.85 ± 3.75</td></tr><tr><td rowspan="2">Fienga et al. (2021a)</td><td rowspan="2">INPOP20a</td><td>PF<sup>**</sup></td><td>-8.8 ± 2.9</td><td>-9.2 ± 6.1</td><td>0.4 ± 9.0</td></tr><tr><td>RW</td><td>-10.3 ± 22.8</td><td>-9.2 ± 6.1</td><td>-0.8 ± 28.4</td></tr></tbody></table>

$$ \frac{\dot{\mu}_P}{\mu_P} = \frac{\dot{G}}{G}, \qquad (137) $$

with $\mu_P$ is the gravitational parameter of bodies other than the Sun.

At each step *t* of the numerical integration of the planetary equations of motion, the following quantities are estimated:

$$ M_{\odot}(t) = M_{\odot}(t_0) + (t - t_0) \times \dot{M}_{\odot}, \qquad (138) $$

$$ M_P(t) = M_P, \qquad (139) $$

$$ G(t) = G(t_0) + (t - t_0) \times \dot{G}, \qquad (140) $$

$$ \mu(t) = G(t) \times M(t), \qquad (141) $$

where $t_0$ is the date of the origin of the planetary ephemeris. Let us note that the gravitational parameter $\mu(t)$ also appears in the computation of the Shapiro delay (see
<table><caption>**Table 6** Values of $\dot{\mu}_{\odot}/\mu_{\odot}$ found in the literature deduced from planetary ephemerides. The Column 3 indicates the type of method used for the estimation: either the full fit (FF), the partial fit (PF) in or the random walk exploration (RW). See text for details. For Fienga et al. (2021a), PF<sup>**</sup> indicates that the oblateness of the Sun is strongly constrained by helio-seismological value. So this solution cannot be seen as a free fit but as a partial fit. Analysis of the MESSENGER mission is labeled MSG in Column 3 using only the Mercury orbit. The $\dot{G}/G$ values are deduced with $\dot{\mu}_{\odot}/\mu_{\odot}$ and $\dot{M}_{\odot}/M_{\odot}$ given in Column 5. The value of $\dot{G}/G$ deduced from Pitjeva and Pitjev (2013) is not the one published by the authors but was obtained in using the same value of $\dot{M}_{\odot}/M_{\odot}$ as for Konopliv et al. (2011)</caption><thead><tr><th>References</th><th>PE</th><th>Method</th><th>μ̇<sub>☉</sub>/μ<sub>☉</sub> × 10<sup>14</sup><br>yr<sup>-1</sup></th><th>Ṁ<sub>☉</sub>/M<sub>☉</sub> × 10<sup>14</sup><br>yr<sup>-1</sup></th><th>GW/G × 10<sup>14</sup><br>yr<sup>-1</sup></th></tr></thead><tbody><tr><td colspan="6"><em>Before AU fixed</em></td></tr><tr><td>Pitjeva (1993)</td><td>EPM1988</td><td>FF</td><td>NA</td><td>NA</td><td>470 ± 470</td></tr><tr><td>Pitjeva (2001)</td><td>EPM2000</td><td>FF</td><td>NA</td><td>NA</td><td>4 ± 8</td></tr><tr><td>Standish (2001)</td><td>DE405</td><td>FF</td><td>NA</td><td>NA</td><td>1 ± 8</td></tr><tr><td>Pitjeva (2005b)</td><td>EPM2004</td><td>FF</td><td>NA</td><td>NA</td><td>-2 ± 5</td></tr><tr><td colspan="6"><em>After AU fixed</em></td></tr><tr><td>Konopliv et al. (2011)</td><td>DE421</td><td>PF</td><td>1 ± 16</td><td>-9.2 ± 6.1</td><td>10.2 ± 22.1*</td></tr><tr><td>Pitjeva and Pitjev (2012)</td><td>EPM2010</td><td>FF</td><td>-5 ± 4</td><td>-6.7 ± 3.1</td><td>1.65 ± 8.77</td></tr><tr><td>Pitjeva and Pitjev (2013)</td><td>EPM2011</td><td>FF</td><td>-6.3 ± 6.4</td><td>-6.7 ± 3.1</td><td>0.4 ± 11.1</td></tr><tr><td rowspan="2">Fienga et al. (2015)</td><td rowspan="2">INPOP15a</td><td>FF</td><td>-5.0 ± 2.9</td><td>-9.2 ± 6.1</td><td>4.2 ± 9.0</td></tr><tr><td>RW</td><td>-4.3 ± 7.4</td><td>-9.2 ± 6.1</td><td>4.9 ± 13.5</td></tr><tr><td>Genova et al. (2018)</td><td>MSG</td><td></td><td>-6.130 ± 1.47</td><td>-10 ± 1</td><td>4 ± 7.5</td></tr><tr><td>Pitjeva et al. (2021)</td><td>EPM2021</td><td>FF</td><td>-10.2 ± 1.4</td><td>-11.05 ± 2.35</td><td>0.85 ± 3.75</td></tr><tr><td rowspan="2">Fienga et al. (2021a)</td><td rowspan="2">INPOP20a</td><td>PF<sup>**</sup></td><td>-8.8 ± 2.9</td><td>-9.2 ± 6.1</td><td>0.4 ± 9.0</td></tr><tr><td>RW</td><td>-10.3 ± 22.8</td><td>-9.2 ± 6.1</td><td>-0.8 ± 28.4</td></tr></tbody></table>

$$ \frac{\dot{\mu}_P}{\mu_P} = \frac{\dot{G}}{G}, \qquad (137) $$

with $\mu_P$ is the gravitational parameter of bodies other than the Sun.

At each step *t* of the numerical integration of the planetary equations of motion, the following quantities are estimated:

$$ M_{\odot}(t) = M_{\odot}(t_0) + (t - t_0) \times \dot{M}_{\odot}, \qquad (138) $$

$$ M_P(t) = M_P, \qquad (139) $$

$$ G(t) = G(t_0) + (t - t_0) \times \dot{G}, \qquad (140) $$

$$ \mu(t) = G(t) \times M(t), \qquad (141) $$

where $t_0$ is the date of the origin of the planetary ephemeris. Let us note that the gravitational parameter $\mu(t)$ also appears in the computation of the Shapiro delay (see
Sect. 3.4.1). In this case, the value of $\mu(t)$, corresponding to the date of the observation, is computed with Eq. (141) and re-introduced in the Shapiro Eq. (94).

$M_\odot(t_0)$ is the mass of the Sun fitted during the construction of the ephemeris, and $G(t_0)$ is the Newtonian gravitation constant as defined by the IAU (Luzum et al. 2011). The effect of the time variation of $G$ being largely induced by the gravitational parameter of the Sun more that the ones of the planets, one then deduces the value of $\dot{G}/G$ by considering its impact of the Sun contribution [Eq. (136)] and a fixed value for the Sun total mass loss, $\dot{M}_\odot/M_\odot$. Pitjeva and Pitjev (2012) have proposed an interval of values for the total solar mass loss of

$$ \frac{\dot{M}_\odot}{M_\odot} = (-0.67 \pm 0.31) \times 10^{-13} \quad (2\sigma) \text{ yr}^{-1}, \qquad (142) $$

This estimation considers the mass loss by the Sun, but also the mass gained by falling materials (comets, asteroids etc...). Pinto et al. (2011) have estimated a mean mass loss from wind emission of charged particles during the 11-year solar cycle:

$$ \frac{\dot{M}_\odot}{M_\odot} = (-0.55 \pm 0.15) \times 10^{-13} \quad (3\sigma) \text{ yr}^{-1}, \qquad (143) $$

instead of Eq. (142). In 2021, a detailed evaluation of Pitjeva et al. (2021) gives

$$ \frac{\dot{M}_\odot}{M_\odot} = (-1.105 \pm 0.235) \times 10^{-13} \quad (3\sigma) \text{ yr}^{-1}, \qquad (144) $$

This value accounts for the solar wind radiation as well as the accumulation of interplanetary dust falling on the Sun together with comets. The authors conclude that the value of $\dot{M}_\odot/M_\odot$ in Pitjeva and Pitjev (2012) is overestimated.

On Table 6, are gathered the values of $\dot{\mu}_\odot/\mu_\odot$ obtained by different authors jointly with PPN parameter estimations as well as the $\dot{G}/G$ deduced using the value of $\dot{M}_\odot/M_\odot$ specified in the same Table. It is important to stress that the correlations between the PPN parameters, $\dot{\mu}_\odot/\mu_\odot$ and the other fitted parameters of the planetary ephemerides being non zero (see e.g. Table 6 from Fienga et al. 2015), the values given in Table 6 have to be considered as part of a global fit, and consequently they are planetary ephemeris dependent. Direct adjustments with all parameters fitted together (Pitjeva and Pitjev 2013; Pitjeva et al. 2021), partial fits (Konopliv et al. 2011; Fienga et al. 2021a) but also random walk exploration algorithms have been used for obtaining constraints on $\dot{\mu}_\odot/\mu_\odot$, given in Table 6. As explained in Sects. 3.5 and 4.1, the partial fit consists to fix one PPN parameter and then fit $\dot{\mu}_\odot/\mu_\odot$ together with the rest of the planetary parameters. The random walk exploration algorithm fixes the PPN parameters and $\dot{\mu}_\odot/\mu_\odot$ according to random values, and performs a regular fit for the rest of the planetary parameters. The obtained ephemerides are then selected according to different statistical criteria (Fienga et al. 2015, 2021a). The latest two approaches (partial fit and random walk exploration algorithms) give larger intervals of possible values than the direct fit of all parameters, as one can see in Table 6.
Sect. 3.4.1). In this case, the value of $\mu(t)$, corresponding to the date of the observation, is computed with Eq. (141) and re-introduced in the Shapiro Eq. (94).

$M_{\odot}(t_0)$ is the mass of the Sun fitted during the construction of the ephemeris, and $G(t_0)$ is the Newtonian gravitation constant as defined by the IAU (Luzum et al. 2011). The effect of the time variation of $G$ being largely induced by the gravitational parameter of the Sun more that the ones of the planets, one then deduces the value of $\dot{G}/G$ by considering its impact of the Sun contribution [Eq. (136)] and a fixed value for the Sun total mass loss, $\dot{M}_{\odot}/M_{\odot}$. Pitjeva and Pitjev (2012) have proposed an interval of values for the total solar mass loss of

$$ \frac{\dot{M}_{\odot}}{M_{\odot}} = (-0.67 \pm 0.31) \times 10^{-13} \quad (2\sigma) \text{ yr}^{-1}, \qquad (142) $$

This estimation considers the mass loss by the Sun, but also the mass gained by falling materials (comets, asteroids etc...). Pinto et al. (2011) have estimated a mean mass loss from wind emission of charged particles during the 11-year solar cycle:

$$ \frac{\dot{M}_{\odot}}{M_{\odot}} = (-0.55 \pm 0.15) \times 10^{-13} \quad (3\sigma) \text{ yr}^{-1}, \qquad (143) $$

instead of Eq. (142). In 2021, a detailed evaluation of Pitjeva et al. (2021) gives

$$ \frac{\dot{M}_{\odot}}{M_{\odot}} = (-1.105 \pm 0.235) \times 10^{-13} \quad (3\sigma) \text{ yr}^{-1}, \qquad (144) $$

This value accounts for the solar wind radiation as well as the accumulation of interplanetary dust falling on the Sun together with comets. The authors conclude that the value of $\dot{M}_{\odot}/M_{\odot}$ in Pitjeva and Pitjev (2012) is overestimated.

On Table 6, are gathered the values of $\dot{\mu}_{\odot}/\mu_{\odot}$ obtained by different authors jointly with PPN parameter estimations as well as the $\dot{G}/G$ deduced using the value of $\dot{M}_{\odot}/M_{\odot}$ specified in the same Table. It is important to stress that the correlations between the PPN parameters, $\dot{\mu}_{\odot}/\mu_{\odot}$ and the other fitted parameters of the planetary ephemerides being non zero (see e.g. Table 6 from Fienga et al. 2015), the values given in Table 6 have to be considered as part of a global fit, and consequently they are planetary ephemeris dependent. Direct adjustments with all parameters fitted together (Pitjeva and Pitjev 2013; Pitjeva et al. 2021), partial fits (Konopliv et al. 2011; Fienga et al. 2021a) but also random walk exploration algorithms have been used for obtaining constraints on $\dot{\mu}_{\odot}/\mu_{\odot}$, given in Table 6. As explained in Sects. 3.5 and 4.1, the partial fit consists to fix one PPN parameter and then fit $\dot{\mu}_{\odot}/\mu_{\odot}$ together with the rest of the planetary parameters. The random walk exploration algorithm fixes the PPN parameters and $\dot{\mu}_{\odot}/\mu_{\odot}$ according to random values, and performs a regular fit for the rest of the planetary parameters. The obtained ephemerides are then selected according to different statistical criteria (Fienga et al. 2015, 2021a). The latest two approaches (partial fit and random walk exploration algorithms) give larger intervals of possible values than the direct fit of all parameters, as one can see in Table 6.