# On the numerical evaluation of the 'exact' Post-Newtonian parameters in Brans-Dicke and Entangled Relativity theories

Thomas Chehab, Olivier Minazzoli

► To cite this version:

Thomas Chehab, Olivier Minazzoli. On the numerical evaluation of the 'exact' Post-Newtonian pa- rameters in Brans-Dicke and Entangled Relativity theories. 2026. hal-05505855

HAL Id: hal-05505855

https://hal.science/hal-05505855v1

Preprint submitted on 11 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On the numerical evaluation of the 'exact' Post-Newtonian parameters in Brans-Dicke and Entangled Relativity theories

Thomas Chehab*

Artemis, Université Côte d'Azur, CNRS, Observatoire Côte d'Azur, BP4229, 06304, Nice Cedex 4, France

Olivier Minazzoli†

Artemis, Université Côte d'Azur, CNRS, Observatoire Côte d'Azur, BP4229, 06304, Nice Cedex 4, France,
Bureau des Affaires Spatiales, 2 rue du Gabian, 98000 Monaco.

In context of Brans-Dicke scalar-tensor theories of gravity, it has recently been obtained that the post-Newtonian parameters should be generalized in the context of strongly gravitating bodies, and that its generalization—the so-called *exact parameters*—actually depends on the pressure and energy density of a considered celestial body. Here we develop two new methods to numerically obtain the *exact parameters* by means of usual Tolman-Oppenheimer-Volkoff computation, and find that the difference with the value of standard post-Newtonian parameters can be more than 80% in some situations. We also provide the connection with the Damour-Esposito Farèse non-pertubative parameter $\alpha_{DEF}$. We then apply the methodology to the case of Entangled Relativity, and derive these exact parameters for the Sun and the Earth, as well as for neutron stars. We argue that current and foreseeable experiments are likely able to constrain the theory under the assumption that $\mathcal{L}_m = -\rho$, where $\rho$ is the total energy density. If $\mathcal{L}_m = T$ instead, as often advocated in the literature, then there is no deviation with respect to General Relativity and the prospects of testing Entangled Relativity become much more remote in time, as only compact objects with extreme electric or magnetic fields could lead to some deviation from General Relativity.

## I. INTRODUCTION

One of the first alternatives to General Relativity was the Brans–Dicke theory of gravity. The main idea behind this theory is the hypothesis that gravity is mediated not only through the metric tensor $g_{\mu\nu}$, but also by a scalar field $\phi$. These theories are known as scalar-tensor theories of gravity [1]. In such frameworks, the scalar field is inversely proportional to the Newtonian gravitational constant $G$. Later, the Bergmann-Wagoner-Nordtvedt theory generalizes the Brans-Dicke theory of gravity using the assumption that the Brans-Dicke coupling parameter $\omega$ can itself depend on the scalar field of the considered scalar-tensor theory [2–4].

Since the first formulation of General Relativity, we have aimed to test the theory experimentally. Firstly, with observations of Mercury’s perihelion, deflection of light by the Sun, and gravitational redshift. These experiment were successful and confirmed to be in agreement with the predictions of General Relativity, giving to it an important credit. Later, the Shapiro delay effect [5] has been studied and nowadays, we impose new constraints using decreasing orbits of binary pulsars [6], gravitational waves [7], direct observations of supermassive Black-hole’s shadow [8], and tests of the validity of the strong equivalence principle with an extreme precision [9]—see [10–13] for reviews. Among all these tests, in the solar system, usual constraints are imposed using post-Newtonian studies. In the weak-field regime,

one can expand the metric in power of $v/c \approx \sqrt{GM/rc^2}$ which leads to coefficients called post-Newtonian parameters. These parameters quantify deviations with respect to General Relativity (for a review of experimental tests of General Relativity, see [14]).

Recently, in the context of Brans–Dicke theory, using the Brans–Dicke class I solution [1], a novel ‘exact’ expression of the post-Newtonian parameters $\beta$ and $\gamma$ has been derived [15] and later in [16], these parameters (plus $\delta$, an higher order parameter) has been derived in Bergmann-Wagoner-Nordvert theory. They has been shown to depend explicitly on the energy density and pressure of the considered body and are exact and non-perturbative quantities, but have so far never been evaluated numerically. Indeed, because these new exact quantities directly depend in the structure of the considered bodies, one needs to integrate the interior solution numerically—in contrast to the case of usual (perturbative) post-Newtonian parameters that can be derived analytically.

In this communication we give a numerical evaluation of these exact (non-perturbative) post-Newtonian parameters in the context of the Brans-Dicke theory and then in Entangled relativity. In order to compute there numerical values we will use the Tolman-Oppenheimer-Volkoff framework following [17–19].

This paper is organized as follows. In Sec. II, we review the notion of exact post-Newtonian parameters in the framework of Bergmann–Wagoner–Nordtvedt scalar-tensor theories and summarize the analytical expressions previously derived for spherical bodies, emphasizing their dependence on the internal structure of the source. In Sec. III, we derive the Tolman–Oppenheimer–Volkoff

<footnote>* chehab.thomas@oca.eu</footnote>

<footnote>† ominazzoli@gmail.com</footnote>
equations in Brans-Dicke theory and explain how these equations can be used to numerically evaluate the exact post-Newtonian parameters. In Sec. IV, we present a detailed numerical study of these parameters in Brans-Dicke gravity, including a cross-check based on external analytical solutions. In addition, we provide the connection with Damour-Esposito Farèse non-pertubative parameter $\alpha_{DEF}$. In Sec. V, we introduce Entangled Relativity, discuss its usual post-Newtonian limit, and derive the exterior solution necessary to compute its exact post-Newtonian parameters. Sec. VI is devoted to the numerical evaluation of these parameters in Entangled Relativity from the exact parameters, for a binary system composed of a pulsar and a white dwarf. Finally, Sec. VIII summarizes our main results.

## II. EXACT POST-NEWTONIAN PARAMETERS IN BERGMANN-WAGONER-NORDTVEDT THEORY

Here we briefly remind the results obtained by [16], in which is derived a generalisation of the exact post-Newtonian parameters in Bergmann-Wagoner-Nordtvedt theory. Its action reads

$$S = \frac{1}{16\pi} \int d^4 x \sqrt{-g} \left[ \Phi R - \frac{\omega(\Phi)}{\Phi} \nabla_{\mu} \Phi \nabla^{\mu} \Phi \right] + S_m[g_{\mu\nu}, \Psi], \quad (1)$$

with $\Phi$ a scalar field, $\omega(\Phi)$ a coupling parameter that depend in the scalar field and $S_m$ the matter action that depend in the metric $g$ and the matter fields $\Psi$. Using the least action principle one can obtain the following field equations

$$G_{\mu\nu} = \frac{8\pi}{\Phi} T_{\mu\nu} + \frac{1}{\Phi} [\nabla_{\mu}\nabla_{\nu} - g_{\mu\nu}\Box]\Phi + \frac{\omega(\Phi)}{\Phi^2} (\nabla_{\mu}\Phi\nabla_{\nu}\Phi - \frac{1}{2}g_{\mu\nu}(\nabla\Phi)^2), \quad (2)$$

$$\Box\Phi = \frac{1}{2\omega(\Phi) + 3} \left[ 8\pi T - \frac{d\omega}{d\Phi} (\nabla\Phi)^2 \right]. \quad (3)$$

### A. Summary of previous results in general scalar-tensor theory

To derive the formulation of the exact post-Newtonian parameters, the author of [16] started in Einstein frame with the Janis-Newman-Winicour metric [20] in isotropic coordinates that they re-expressed in the Jordan frame thanks to a conformal transformation. This latter reads

$$d\varphi = \sqrt{3 + 2\omega(\Phi)} d \ln \Phi. \quad (4)$$

While the metric solution reads

$$ds^2 = -Adt^2 + Bdr^2 + Cd\Omega^2, \quad (5)$$

Where the function A, B and C are given by

$$\left\{ \begin{array}{l} A = \frac{1}{\Phi} \left( \frac{r-k}{r+k} \right)^{2d}, \\ B = \frac{1}{\Phi} \left( 1 + \frac{k}{r} \right)^4 \left( \frac{r-k}{r+k} \right)^{2-2d}, \\ C = Br^2. \end{array} \right. \quad (6)$$

In which $k$ is linked to the mass by $k = \frac{2m}{d}$ and $d$ is the scalar charge. They assumed a static and spherically symmetric stress energy tensor defined by

$$(T_a^b) = \operatorname{diag}(-\rho, p_\parallel, p_\perp, p_\perp). \quad (7)$$

With $\rho$ the energy density and $p_\parallel$, $p_\perp$ are respectively the parallel and perpendicular pressure. Both are functions of $r$. By focusing on the (00) component of Eq. (2) and the scalar field equation Eq. (3) one can derive the value of $d$ and $k$ as

$$\left\{ \begin{array}{l} d = \frac{E^* + P^*}{\sqrt{(E^* + P^*)^2 + (E_\omega^* - P_\omega^*)^2}}, \\ k = \frac{1}{2} \sqrt{(E^* + P^*)^2 + (E_\omega^* - P_\omega^*)^2}, \end{array} \right. \quad (8)$$

with

$$\left\{ \begin{array}{l} E^* = \int_{(\text{star})} (-T_0^0) dV, \\ P^* = \int_{(\text{star})} T_i^i dV, \\ E_\omega^* = \int_{(\text{star})} \frac{(-T_0^0)}{\sqrt{3+2\omega(\Phi)}} dV, \\ P_\omega^* = \int_{(\text{star})} \frac{T_i^i}{\sqrt{3+2\omega(\Phi)}} dV, \end{array} \right. \quad (9)$$

where $E^*$, and $P^*$ are respectively the total integrated energy and total integrated pressure while $E_\omega^*$, and $P_\omega^*$ are the total integrated energy and pressure divided by a function of $\omega(\Phi)$. Here $dV = 4\pi\sqrt{ABC}dr$. We see with Eq. (8) that the scalar charge $d$ explicitly depend on the structure of the considered body.

By Taylor expanding A and B one can obtain

$$\left\{ \begin{array}{l} A = 1 + x \left[ \frac{d}{dx}(A) \right]_0 + \frac{1}{2} x^2 \left[ \frac{d^2}{dx^2}(A) \right]_0 + O(x^3), \\ B = 1 + x \left[ \frac{d}{dx}(B) \right]_0 + \frac{1}{2} x^2 \left[ \frac{d^2}{dx^2}(B) \right]_0 + O(x^3), \end{array} \right. \quad (10)$$

with $x$ given by $x = \frac{k}{r}$. The subscript 0 indicates that the value is taken at $x = 0$ (ie : $r \to \infty$). By solving the second and third term of A and B one can compute the total metric expansion and thus identify, by comparing with usual post-Newtonian definition, the mass and the exact post-Newtonian parameters as
$$m = 2 \left( d + \sqrt{\frac{1-d^2}{3+2\omega_0}} \right) k, \quad (11a)$$

$$\beta_{\text{exact}} = 1 + \left( \frac{\left(\frac{d\omega}{d\Phi}\right)_0}{(3+2\omega_0)^2} \right) \frac{1-d^2}{\left( d + \sqrt{\frac{1-d^2}{3+2\omega_0}} \right)^2}, \quad (11b)$$

$$\gamma_{\text{exact}} = \frac{d - \sqrt{\frac{1-d^2}{3+2\omega_0}}}{d + \sqrt{\frac{1-d^2}{3+2\omega_0}}}, \quad (11c)$$

$$\delta_{\text{exact}} = \frac{4}{3} (\beta_{\text{exact}} + \gamma_{\text{exact}}^2 - 1) - \quad (11d)$$

$$\frac{1}{3} \left( d + \sqrt{\frac{1-d^2}{3+2\omega_0}} \right)^{-2},$$

where $\omega_0$ is the background value of $\omega(\Phi)$. As a reminder, $d$ is a parameter of the Janis-Newman-Winicour metric (Eq.(5)) and is given by Eq.(8). It explicitly depends on the structure of the star. Hence the post-Newtonian parameters $\beta_{\text{exact}}$, $\gamma_{\text{exact}}$, $\delta_{\text{exact}}$ are formulated with an explicit dependence in the structure of the considered body. This results has been obtained using isotropic Janis-Newman-Winicour metric. Using this result one can replace $d$ with Eq.(8) an re-arrange the equation to obtain

$$m = \left( 1 + \frac{\Xi}{\sqrt{3 + 2\omega_0}} \right) (E^* + P^*), \quad (12a)$$

$$\frac{\sigma}{m} = \frac{\Xi}{\sqrt{3 + 2\omega_0 + \Xi}}, \quad (12b)$$

$$\beta_{\text{exact}} = 1 + \frac{\left(\frac{d\omega}{d\Phi}\right)_0}{3 + 2\omega_0} \left(\frac{\sigma}{m}\right)^2, \quad (12c)$$

$$\gamma_{\text{exact}} = \frac{\sqrt{3 + 2\omega_0} - \Xi}{\sqrt{3 + 2\omega_0} + \Xi}, \quad (12d)$$

$$\delta_{\text{exact}} = \frac{4}{3} \left[ \beta + \gamma^2 - 1 - \frac{(3 + 2\omega_0)(1 + \Xi^2)}{4 \left(\sqrt{3 + 2\omega_0} + \Xi\right)^2} \right] \quad (12e)$$

in which $\Xi = \frac{1}{\sqrt{3+2\omega_0}} \frac{1-\theta}{1+\theta}$ with $\theta = \frac{P^*}{E^*}$. Detailed calculation are presented in [16].

## B. Invariance under radial coordinate transformation

An important point lies in the fact that the integrals presented in equations (9) are invariant under a radial coordinate change. Indeed if one practices such a change $(t, r, \theta, \phi) \rightarrow (t' = t, r'(r), \theta' = \theta, \phi' = \phi)$ in the end, the volume element gives: $\sqrt{A'B'C'dr'} = \sqrt{AB \left(\frac{dr}{dr'}\right)^2 Cdr'}$. And hence $dV = dV'\text{sign}(\frac{dr}{dr'})$. Then, the sign of $\frac{dr}{dr'}$ can be reabsorbed into the permutation of lower and upper bound in the integral. Lastly, $T_0^0$ and $T_i^i$ are scalars.

Hence these integrals are invariant through a radial coordinate change.

This is important because the Tolman-Oppeheimer-Volkoff metric does not use the same radial coordinates than the one used by [16]. Hence, this invariance is a necessary property in order to be able to inject the Tolman-Oppeheimer-Volkoff numerical evaluations of the various quantities in the analytical expressions provided in [16]—summarized in Sec. II A.

## C. Difference between exact and usual post-Newtonian parameters

It is worth emphasizing that the exact post-Newtonian parameters in equations (11-12) depend on the internal structure of each bodies, whereas the usual post-Newtonian parameters in the literature, see e.g. [10, 12], only depend on the specific theory considered.

The difference between the two set of post-Newtonian parameters is that the usual post-Newtonian parameters are derived in the weak-field limit of a given theory—in particular, assuming that $\rho \gg P$—whereas the exact post-Newtonian parameters are non-perturbative quantities that are also valid for ultra-relativistic compact bodies—for which $P \sim \rho$.

The issue with exact parameters in equations (11-12) is that one needs to know the internal structure of the star—through $\rho$, $p_\parallel$ and $p_\perp$—as well as the internal metric—through $A$ and $B$—in order to evaluate it. In what follows, we will show how to obtain it numerically by solving the Tolman-Oppeheimer-Volkoff equations.

## III. TOLMAN-OPPENHEIMER-VOLKOFF EQUATIONS IN BRANS-DICKE THEORY - REWRITING THE EXACT POST-NEWTONIAN PARAMETERS

In the Bergmann-Wagoner-Nordtvedt theory, the parameter $\omega$ is a function of $\phi$. As a consequence, the Bergmann-Wagoner-Nordtvedt theory actually corresponds to an infinite class of theories—one for each function $\omega(\Phi)$ considered—and one needs to restrict oneself to a specific theory in order to derive the corresponding phenomenology. We will focus in the particular case $\omega(\phi) = \omega$ which constitute the case of the Brans-Dicke theory. Current constraints on the Strong Equivalence Principle with a pulsar in a triple star system gives $\omega > 140 000$ [21].

We will firstly derive the modified Tolman-Oppenheimer-Volkoff equations of the Brans-Dicke theory that we will solve numerically in order to evaluate the exact post-Newtonian parameters. Then, with numerical integration we will directly evaluate these exact parameters.
A. Tolman-Oppenheimer-Volkoff equations of the Brans-Dicke theory

To derive the Tolman-Oppenheimer-Volkoff equations in Brans-Dicke theory, we will follow the same method as in [17]. The Brans-Dicke theory [1] is defined by the action

$$S_{BD} = \int d_g^4 x \left[ \frac{1}{2\tilde{\kappa}} \left( \phi R - \frac{\omega}{\phi} \partial_\sigma \phi \partial^\sigma \phi \right) + \mathcal{L}_m \right]. \quad (13)$$

Using the least action principle one can obtain the following field equations

$$G_{\mu\nu} = \frac{\tilde{\kappa}}{\phi} T_{\mu\nu} + \frac{\omega}{\phi^2} \left( \partial_{\mu}\phi\partial_{\nu}\phi - \frac{1}{2}g_{\mu\nu}\partial_{\sigma}\partial^{\sigma}\phi \right) + \frac{1}{\phi} \left( \nabla_{\mu}\nabla_{\nu}\phi - g_{\mu\nu}\Box\phi \right), \quad (14)$$

$$\Box\phi = \frac{\tilde{\kappa}}{3 + 2\omega} T, \quad (15)$$

with $\tilde{\kappa} = \frac{8\pi G}{c^4}$, with G the Newton constant. We now set:

$$D_{\beta}^{\alpha} = \frac{1}{\phi} \left( \nabla^{\alpha} \nabla_{\beta} - \delta_{\beta}^{\alpha} \Box \right) \phi, \quad (16)$$

and

$$H_{\beta}^{\alpha} = \frac{\omega}{\phi^2} \left( \partial^{\alpha} \phi \partial_{\beta} \phi - \frac{1}{2} \delta_{\beta}^{\alpha} \partial_{\sigma} \partial^{\sigma} \phi \right). \quad (17)$$

Using the (00) component of the field equation (14), we obtain an equation over $g_{rr} \equiv b$.

$$\frac{\dot{b}}{b} = rb \left( \frac{\tilde{\kappa}\rho}{\phi} - H_0^0 - D_0^0 \right) + \frac{1-b}{r}. \quad (18)$$

While one finds:

$$D_0^0 = \frac{\dot{a}}{2ab} \frac{\dot{\phi}}{\phi} + \frac{\tilde{\kappa}}{3+2\omega} \frac{T}{\phi}, \quad (19)$$

and

$$H_0^0 = -\frac{\omega \dot{\phi}^2}{2b \phi^2}, \quad (20)$$

in which $a \equiv g_{tt}$. With a direct integration of Eq. (18), one obtains

$$b = \frac{1}{1 - \frac{m\tilde{\kappa}}{4\pi r}}, \quad (21)$$

setting

$$m = \frac{4\pi}{\tilde{\kappa}} \int_{0}^{r} dr' r'^2 x \left( \frac{\tilde{\kappa}\rho}{\phi} - H_{0}^{0} - D_{0}^{0} \right). \quad (22)$$

Thanks to which we can derive a differential equation for the mass

$$\dot{m} = \left( \frac{\tilde{\kappa}\rho}{\phi} - H_0^0 - D_0^0 \right) \frac{4\pi r^2}{\tilde{\kappa}}. \quad (23)$$

Doing the same work with the (11) composant of the field equation (14), leads to differential equation for a

$$\begin{aligned} \frac{\dot{a}}{a} = \frac{b}{r} & \left( \frac{r^2 \tilde{\kappa} P}{\phi} + \frac{\dot{\phi}^2}{\phi^2} \frac{r^2 \omega}{2b} + 1 \right) \left( 1 + \frac{r}{2} \frac{\dot{\phi}}{\phi} \right)^{-1} \\ & - \frac{b}{r} \left( \frac{1}{b} + \frac{2r}{b} \frac{\dot{\phi}}{\phi} \right) \left( 1 + \frac{r}{2} \frac{\dot{\phi}}{\phi} \right)^{-1} . \end{aligned} \quad (24)$$

Concerning the pressure, using the conservation equation of the theory $\nabla_\sigma T^{\alpha\sigma} = 0$, one can obtain a differential equation representing its behaviour

$$\dot{P} = -(\rho + P) \frac{\dot{a}}{2a}. \quad (25)$$

Lastly, with the scalar field equation (15), one obtains a differential equation of the scalar field

$$\ddot{\phi} = -\frac{\dot{\phi}}{2} \left( \frac{\dot{a}}{a} - \frac{\dot{b}}{b} + \frac{4}{r} \right) + b \Box \phi. \quad (26)$$

a. Summary of equations: Up to now, every necessary quantity has been expressed as a differential equation that can be integrated numerically to model compact objects in Brans-Dicke theory.

$$\left\{ \begin{aligned} \dot{P} &= -(\rho + P) \frac{\dot{a}}{2a}, \\ \dot{m} &= \left( \frac{\tilde{\kappa}\rho}{\phi} - H_0^0 - D_0^0 \right) \frac{4\pi r^2}{\tilde{\kappa}}, \\ \dot{\phi} &= \psi, \\ \dot{\psi} &= -\frac{\psi}{2} \left( \frac{\dot{a}}{a} - \frac{\dot{b}}{b} + \frac{4}{r} \right) + b \Box \phi. \end{aligned} \right. \quad (27)$$

B. The equation of state

The goal of this study is to numerically evaluate exact post-Newtonian parameters. Hence our goal is only to obtain their broad behaviour. This is why as a first exploratory model, we decided to use a simple polytropic equation of state defined by
$$P = K \rho^{\gamma}, \qquad (28)$$

with $\gamma = 5/3$, is the adiabatic index and $K = 1.475 \times 10^{-3}$ (fm$^3$/MeV)$^{2/3}$ is simply a constant [17, 22]. Moreover, we also do not take into account a potential effect of a magnetic field or a rotation, for the same reason as explain above.

### C. Solving the Boundary Value Problem

Solving these equations is a Boundary Value Problem. Indeed, one wants to find the initial condition at the center of the star that will lead to the desired normalization at infinity for any value of the central density. To resolve this problem we used a normalization technique [17]. Indeed one can remark that if $(\phi, \psi, m, P)(r)$ are solution of equation (27), then another solution is given by $(\alpha\phi, \alpha^{5/4}\psi, \alpha^{-1/4}m, \alpha^{3/2}P)(\alpha^{-1/4}r)$ for any $\alpha$. Then if a simulation starts with an initial value $\phi_0$ that leads to a given star radius $R_*$ with an asymptotical scalar field value $\phi_\infty$, then the radius with the appropriate unit is given by $R_* = \phi_\infty^{1/4}R_*$, when the property has been used with $\alpha = \phi_\infty$, such that the final normalization of $\phi$ at infinity is one.

### D. Exact post-Newtonian parameter in Brans-Dicke gravity

As explain in sec. IIB, the equations for the exact post-Newtonian parameters given in Eq. (12) are invariant under radial coordinate change. Hence even if they have been derived in [16] using Janis-Newman-Winicour metric in isotropic coordinate, the result is the same in a Tolman-Oppenheimer-Volkoff spherical metric.

In the framework of Brans dicke theory, as $\omega = \text{cst}$ one can rewrite Eq. (12) as

$$\beta_{\text{exact}} = 1, \qquad (29a)$$

$$\gamma_{\text{exact}} = \frac{1 + \omega + \theta(2 + \omega)}{2 + \omega + \theta(1 + \omega)}, \qquad (29b)$$

$$\delta_{\text{exact}} = \frac{4}{3} \left[ \gamma_{\text{exact}}^2 - \frac{(3 + 2\omega)(1 + \Xi^2)}{4 \left(\sqrt{3} + 2\omega + \Xi\right)^2} \right]. \qquad (29c)$$

With $\Xi = \frac{1}{\sqrt{3+2\omega}} \left(\frac{1-\theta}{1+\theta}\right)$ in Brans-Dicke theory. As a reminder $\theta = \frac{P^*}{E^*}$, with

$$X^* = \int_0^{\rho^*} 4\pi X \sqrt{ab^3} \rho^2 d\rho, \qquad (30)$$

where $\rho$ is the radial coordinate in Tolman-Oppenheimer-Volkoff framework. Then, numerical evaluations of $a, b, \rho$ and $P$ allows one to evaluate Eq. (30) numerically, from which one infers the values of the exact post-Newtonian parameters in Eq. (29).

## IV. NUMERICAL EVALUATION OF EXACT PARAMETERS IN BRANS-DICKE THEORY

### A. Numerical estimation of $(\gamma_{\text{exact}}, \beta_{\text{exact}}, \delta_{\text{exact}})$ for neutron stars

We can now evaluate exact post-Newtonian parameters numerically. Indeed using numerical integration, we solve the Tolman-Oppenheimer-Volkoff equations in the framework of the Brans-Dicke theory.

In figures 1 and 2 we have plotted a shade of colour that represents $1 - \gamma_{\text{exact}}$ and $\frac{\gamma_{\text{exact}} - \gamma_{PN}}{1 - \gamma_{PN}} \%$, respectively, to probe there respective deviations for various values of the parameter $\omega$ and the central energy density $\rho_c$. In figures 3 and 4, the shade of colour represents respectively $1 - \delta_{\text{exact}}$ and $\frac{\delta_{\text{exact}} - \delta_{PN}}{1 - \delta_{PN}} \%$. Here $\gamma_{PN}$ and $\delta_{PN}$ are the standard (perturbative) formulation of the post-Newtonian parameters $\gamma$ and $\delta$ and are given by $\gamma_{PN} = \frac{1+\omega}{2+\omega}$ and $\delta_{PN} = \frac{4}{3}(\gamma_{PN}^2 + \beta_{PN}) - \frac{\gamma}{6} - \frac{3}{2}$ [23, 24].

One can notice in all figures, that for a same value of $\omega$, the resulting colors values explicitly varies with respect to the density of the star, which is expected as the exact post-Newtonian parameters depend in the structure of the object studied.

<figure><img src="image_2.png" /><figcaption>FIG. 1: Plot of $1 - \gamma_{\text{exact}}$ for 150 different values of $\omega$ and of central densities. The values of $\omega$ ranges from $10^{-1}$ to $10^6$ meanwhile the values of central density ranges from 100 MeV/fm$^3$ to 1570 MeV/fm$^3$. This limit comes from the fact that we want the speed of sound in the neutron star to be less than the conservative limit considered of $c/\sqrt{3}$ [25]. The deviation of $\gamma_{\text{exact}}$ from unity is represented in colour shades.</figcaption></figure>

In Fig. 1 we see that the deviation from unity of $\gamma_{\text{exact}}$ decreases with increasing $\omega$ and also with increasing energy density. This is expected for increasing $\omega$, but it is somewhat counter intuitive for increasing central energy densities, given that $\Phi \propto T$. But let us recall that the scalar field equations reads
<figure><img src="image_3.png" /><figcaption>FIG. 2: Plot of the relative difference between the exact and standard post-Newtonian parameters $\gamma$ for 150 different values of $\omega$ and central densities. The values of $\omega$ ranges from $10^{-1}$ to $10^6$ meanwhile the values of central density ranges from 100 MeV/fm$^3$ to 1570 MeV/fm$^3$. This limit comes from the fact that we want the speed of sound in the neutron star to be less that the conservative limit considered of $c/\sqrt{3}$ [25]. In color shades is represented the relative deviation in percent of both $\gamma$ definitions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3: Plot of $1 - \delta_{\text{exact}}$ for 150 different values of $\omega$ and central densities. The values of $\omega$ ranges from $10^{-1}$ to $10^6$ meanwhile the values of central density ranges from 100 MeV/fm$^3$ to 1570 MeV/fm$^3$. This limit comes from the fact that we want the speed of sound in the neutron star to be less that the conservative limit considered of $c/\sqrt{3}$ [25]. In color shades is represented $\delta_{\text{exact}}$.</figcaption></figure>

$$ \Box\phi = \frac{8\pi}{3+2\omega} (-\rho + p_{\parallel} + 2p_{\perp}), \quad (31) $$

such that

$$ \Box\phi \approx \frac{8\pi\rho}{3+2\omega} \left(-1+\frac{3p}{\rho}\right). \quad (32) $$

Therefore, because $P \to \rho/3$ for the most compact objects for which the speed of sound is close to $c/\sqrt{3}$,<fnref n="1" /> the term in the parenthesis becomes closer to zero with an increasing central density. This explains why, somewhat paradoxically, fewer deviation from General Relativity is found for the most compact objects. Overall, the deviation from unity is found to be up to 44.1%.

In Fig. 2, one can see that, as the density goes up, more and more difference appears between the exact and standard post-Newtonian parameters $\gamma$. Indeed, for a given value of $\omega$, the value of $\gamma_{BD}$ stay the same, whereas the one of $\gamma_{\text{exact}}$ approaches the General Relativity value—because of the cancellation in Eq. (32) in the $P \to \rho/3$ limit.

Hence, the relative variation increases with density. The relative variation is found to reach more than 44.7%.

In Fig. 3, we see the same behavior for $1 - \delta_{\text{exact}}$, as the scalar field is less and less sourced with increasing $\omega$ and increasing $\rho$. The deviation from unity is found to be up to 83.7%.

Lastly, in Fig. 4, one can see that, as the density goes up, more and more difference appears between exact and standard post-Newtonian parameters $\delta$. For a given value of $\omega$, the value of $\delta_{BD}$ stay the same whereas the one of $\delta_{\text{exact}}$ approaches the General Relativity value—because of the cancellation in Eq. (32) in the $P \to \rho/3$ limit. Hence the relative variation increases with density. The relative variation is found to reach 46.3%.

The code that generates all the figures presented in this manuscript is freely available on github [26].

## B. Verification by comparing to an external numerical solution

In order to verify our results we will use an exact analytical solution derived in [18]. The method is to start from a general Janis-Newman-Winicour metric and conformally transform it to the Jordan frame. Thanks to what we can recover the Tolman-Oppenheimer-Volkoff metric. In what follows we explain how this works.

<footnote n="1">¹ Let us recall that the speed of sound $c_s$ is defined locally as $(c_s/c)^2 = \partial P/\partial\rho$, such that $c_s = c/\sqrt{3}$ when $P = \rho/3$. Let us also recall that, following [25], we limited our integrations to cases that are such that $c_s = c/\sqrt{3}$.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4: Plot of the relative difference between the exact and standard post-Newtonian parameters $\delta$ for 150 different values of $\omega$ and central densities. The values of $\omega$ ranges from $10^{-1}$ to $10^6$ meanwhile the values of central density ranges from 100 MeV/fm$^3$ to 1570 MeV/fm$^3$. This limit comes from the fact that we want the speed of sound in the neutron star to be less that the conservative limit considered of $c/\sqrt{3}$ [25]. In color shades is represented the relative deviation in percent of both $\delta$ definitions.</figcaption></figure>

The Brans-Dicke action in the Einstein frame reads

$$S_{BD} = \frac{1}{2\tilde{\kappa}} \int d_{\tilde{g}}^4 x \left[ \tilde{R} - 2 (\partial\varphi)^2 + \tilde{\mathcal{L}}_m (\varphi, g_{\mu\nu}; \Psi) \right] \quad (33)$$

with the field definition<fnref n="2" />

$$d\varphi = \pm \frac{\sqrt{3+2\omega}}{2} d \ln \phi. \quad (34)$$

A class of vacuum spherical solutions is given by the Janis-Newman-Winicour metric:

$$ds^2 = - \left(1 - \frac{k}{r}\right)^d dt^2 + \left(1 - \frac{k}{r}\right)^{-d} dr^2 + r^2 \left(1 - \frac{k}{r}\right)^{1-d} d\Omega^2, \quad (35)$$

and the scalar field is given by

$$e^{\varphi} = \left(1 - \frac{k}{r}\right)^{\frac{1-d}{2\alpha}}, \quad (36)$$

with $\alpha^2 = \frac{1-d^2}{1+d^2}$. Here, $d$ ranges between $[-1; 0[\cap]0; 1]$ and $d\Omega^2 = d\theta^2 + \sin^2(\theta)d\phi^2$ is the 2-sphere element. This

writing of the Janis-Newman-Winicour solution can already be found in [23].

For later convenience, it is useful to rewrite the solution using

$$d = \frac{1 - \alpha^2}{1 + \alpha^2}, \quad (37)$$

so that, one obtains

$$ds^2 = - \left(1 - \frac{k}{r}\right)^{\frac{1-\alpha^2}{1+\alpha^2}} dt^2 + \left(1 - \frac{k}{r}\right)^{-\frac{1-\alpha^2}{1+\alpha^2}} dr^2 + r^2 \left(1 - \frac{k}{r}\right)^{\frac{2\alpha^2}{1+\alpha^2}} d\Omega^2. \quad (38)$$

With $\alpha \in \mathbb{R} - \{1\}$, and

$$e^{\varphi} = \left(1 - \frac{k}{r}\right)^{\frac{\alpha}{1+\alpha^2}}. \quad (39)$$

We can now perform the inverse conformal transformation to set this Janis-Newman-Winicour metric into the Jordan frame to recover Tolman-Oppenheimer-Volkoff metric to study this solution numerically. The inverse transformation is given by

$$g_{\alpha\beta} = e^{\mp\xi\varphi} \tilde{g}_{\alpha\beta}, \quad (40)$$

with $\phi = e^{\mp\xi\varphi}$ in which $\xi$ is the conformal factor given by $\xi = \frac{2}{\sqrt{3+2\omega}}$ in the context of the Brans-Dicke theory. The resulting metric reads

$$ds^2 = -adt^2 + bd\rho^2 + \rho^2 d\Omega^2. \quad (41)$$

With

$$a = \left(1 - \frac{k}{r}\right)^{\frac{1-\alpha^2 \mp \alpha \xi}{1+\alpha^2}}, \quad (42a)$$

$$b = \left(\frac{dr}{d\rho}\right)^2 \left(1 - \frac{k}{r}\right)^{\frac{\alpha^2 - 1 \mp \alpha \xi}{1+\alpha^2}}, \quad (42b)$$

$$\rho = r \left(1 - \frac{k}{r}\right)^{\frac{2\alpha^2 \mp \alpha \xi}{2(1+\alpha^2)}}, \quad (42c)$$

$$\phi = \left(1 - \frac{k}{r}\right)^{\pm \frac{\alpha \xi}{1+\alpha^2}}. \quad (42d)$$

Hence, using conformal transformation of the Brans-Dicke theory, starting from the Janis-Newman-Winicour metric, we recover Tolman-Oppenheimer-Volkoff metric. We have thus two ways to compute the exact parameters.

<footnote n="2">² This normalization is different from the one used by [16], and was chosen to be consistent with other previous works [27–30].</footnote>
a. **Computing exact parameters using numerical evaluation.** This is the method described in Sec. III D. In this case we explicitly use numerical evaluation of exact parameters<fnref n="3" /> (equations (29)) by integrating numerically Eq. (30).

b. **Computing exact parameters using external solution.** Another way comes from the evaluation of the parameter $d = \frac{1-\alpha^2}{1+\alpha^2}$ of the Janis-Newman-Winicour solution, from the external part of numerical solutions. Indeed using Tolman-Oppenheimer-Volkoff solution given by Eq. (41), one obtain equation that governs the scalar field, in the vacuum limit

$$ \frac{\ddot{\phi}}{\dot{\phi}} = \frac{1}{2} \left[ \frac{\dot{b}}{b} - \frac{\dot{a}}{a} - \frac{4}{\rho} \right]. \qquad (43) $$

Then, by integration, one obtains

$$ \dot{\phi} = \frac{C}{\rho^2} \sqrt{\frac{b}{a}}, \qquad (44) $$

with $C$ a constant of integration. Finally using Eq. (39) and Eq. (42), by computing the limit when $\rho \to \infty$, one obtains

$$ \alpha = \frac{-\frac{2}{\sqrt{3}}(1+C) - \sqrt{\Delta}}{2}, \qquad (45) $$

with $\Delta$ given by $\Delta = \left(\frac{2}{\sqrt{3}}(1+C)\right)^2 + 4$.

Thanks to which we can evaluate exact post-Newtonian parameters, in Brans-Dicke theory (i.e., by switching from $\alpha$ to $d$ with Eq. (37) and applying $\omega = \text{cst}$ to equations (11)) given by

$$ \beta_{\text{exact}} = 1, \qquad (46a) $$

$$ \gamma_{\text{exact}} = \frac{d - \sqrt{\frac{1-d^2}{3+2\omega_0}}}{d + \sqrt{\frac{1-d^2}{3+2\omega_0}}}, \qquad (46b) $$

$$ \delta_{\text{exact}} = \frac{4}{3}\gamma_{\text{exact}}^2 - \qquad (46c) $$

$$ \frac{1}{3} \left( d + \sqrt{\frac{1-d^2}{3+2\omega_0}} \right)^{-2} . $$

The method that deduces the Janis-Newman-Winicour parameter $\alpha$ from the external part of numerical Tolman-Oppenheimer-Volkoff solutions was originally presented in [18].

c. **Cross-checking results.** Hence we have two methods to compute the exact parameters, such that we can cross-check their results. One that fits the Janis-Newman-Winicour parameter $\alpha$ (or $d$) to the external part of numerical solutions, and the other that numerically computes analytical expressions that depends on $a$, $b$, $\rho$ and $P$ originally presented in [16]—see Sec. III D.

We have verified that the highest relative deviation in absolute value, computed by our code between both $\gamma_{\text{exact}}$ estimations are below 0.07% and at most of 0.3% for $\delta_{\text{exact}}$.

This notably is an independent test of the validity of the analytical results presented in [16].

### C. Discussion on Bergmann-Wagoner-Nordvert theories

General Bergmann-Wagoner-Nordvert theories are known to have non-pertubative strong-field effects that enhance the scalar charge of compact objects [31]. This process is known as the ‘spontaneous scalarization’ of compact objects [32]. It has notably been used in order to constrain scalar-tensor theories, notably with white dwarf-neutron star systems, through the fact that a strong dipolar emission of gravitational waves is expected from scalarized compact objects, whereas no sign of dipolar emission is found in the observational data to date [33].

Therefore, since the exact post-Newtonian parameters are non-pertubative quantities, one can expect spontaneous scalarizations to induce a significant deviation of the exact post-Newtonian parameters from the standard post-Newtonian parameters. This should then be compared to the body-dependent parameters $\alpha_A$ and $\beta_A$ in the Damour-Esposito-Farèse parametrization of Scalar-Tensor theories of gravity [32, 34], which are also sensitive to non-pertubative effects—unlike standard post-Newtonian parameters. This is the topic of the following sub-section.

### D. Linking Damour-Esposito Farèse non-pertubative parameter to $\gamma_{\text{exact}}$

One can expect a direct link between the scalar charge of compact object and exact post-Newtonian parameters. Non-perturbative strong field effects have notably been studied in binary system by [35]. Starting from a scalar-tensor theory in Jordan frame, an arbitrary coupling function $A(\varphi)$ arise when changing to the Einstein frame as $\tilde{g}_{\mu\nu} = A^2(\varphi)g_{\mu\nu}$. Following [35], in the far scalar field region we have

$$ \varphi = \varphi_0 - \frac{m_A \alpha_{DEF}}{r} + \mathcal{O}(r^{-2}), \qquad (47) $$

<footnote n="3">³ Let us note that while we use different radial coordinate systems in Sec. II and Sec. IV, the numerical evaluation of all the relevant quantities is invariant under this change, as emphasized in Sec. II B.</footnote>
in which $\varphi$ is the scalar field and $\varphi_0$ its value at infinity. $m_X$ represent the mass of the body X and $\alpha_{DEF}(\varphi) = \frac{\partial \ln A(\varphi)}{\partial \varphi}$. The subscript DEF stands for Damour and Esposito-Farèse. This should be compared to the expansion of the scalar field equation (39) of the Janis-Newman-winicour solution, in the far field region. This one reads

$$ \varphi = \varphi_0 - \frac{\alpha_{JNW}}{1 - \alpha_{JNW}^2} \frac{2m_A}{r} + \mathcal{O}(r^{-2}). \quad (48) $$

The subscript JNW stands for Janis-Newman-Winicour. By identification, we have

$$ \alpha_{DEF} = \frac{2\alpha_{JNW}}{1 - \alpha_{JNW}^2}. \quad (49) $$

Hence we can derive a link between the scalar charge of Damour-Esposito Farèse and the one of Janis-Newman-Winicour. Using Eq. (11c) and keeping in mind Eq. (37), one finally obtains a formulation depending on the exact Post-Newtonian parameter, for the Brans-Dicke theory,

$$ \alpha_{DEF} = \sqrt{3 + 2\omega_0} \left( \frac{1 - \gamma_{\text{exact}}}{1 + \gamma_{\text{exact}}} \right), \quad (50) $$

or, alternatively,

$$ \gamma_{\text{exact}} = \frac{\sqrt{3 + 2\omega_0} - \alpha_{DEF}}{\sqrt{3 + 2\omega_0} + \alpha_{DEF}}. \quad (51) $$

## V. EXACT PARAMETERS IN ENTANGLED RELATIVITY

The aim of this section is to provide a brief overview of the theory of Entangled Relativity and then to examine the standard post-Newtonian expansion of the theory and derive an exterior solution of it, which will be useful in the next section for the numerical evaluation of $\gamma_{\text{exact}}$ and $\delta_{\text{exact}}$ within the solar system and neutron stars.

### A. The theory

One of the motivation of this theory comes from the fact that by a simple modification of the Einstein-Hilbert action one can immediately satisfy Mach's principle<fnref n="4" /> while being more economical in terms of fundamental

constant than General Relativity and still recovering the phenomenology of General Relativity in some limit. Mach's principle—or principle of relativity of inertia—states that movement, and thus inertia, cannot be defined with respect to nothing. In relativistic theories, inertia is defined through the metric tensor. But General Relativity admits vacuum solutions, thereby automatically violating Mach's principle.

In Entangled Relativity, instead of the standard additional coupling of General Relativity between matter and spacetime in the action, they are coupled non-linearly—see Eq. (52). Because of that, the theory cannot even be defined without matter fields. Moreover, the field equations are such that vacuum solutions are ill-defined, and therefore prohibited—although it has been shown that all the vacuum solutions of General Relativity are limits of solutions of Entangled Relativity when the density of matter fields tends toward zero [40, 41].

To see this more clearly, let us define the theory from its action [19, 42–44]

$$ S_{ER} \propto \int d^4 x \sqrt{-g} \frac{\mathcal{L}_m^2}{R}, \quad (52) $$

where R is the usual Ricci scalar that is constructed upon the metric tensor $g$, $d_g^4 x := \sqrt{-|g|} d^4 x$ is the spacetime volume element, with $|g|$ the metric $g$ determinant, and $\mathcal{L}_m$ is the Lagrangian density of matter. As one can see in Eq. (52), if $\mathcal{L}_m(f, g) = \emptyset$ in the action, the theory cannot be defined at all.

It has been shown in [42, 44, 45] that there is a one to one correspondence at the classical level between the action of Entangled Relativity and an Einstein-dilaton theory (provided that $\mathcal{L}_m \neq \emptyset$) with the following action

$$ S_{Ed} \propto \int d^4 x \sqrt{-g} \left( \frac{\phi R}{2\tilde{\kappa}} + \sqrt{\phi} \mathcal{L}_m \right). \quad (53) $$

Here, $\tilde{\kappa}$ is a positive effective coupling constant between matter and geometry, with the dimension of $\kappa$ from General Relativity and $\phi$ is a dimensionless scalar field. The classical equivalence between the original $f(R, \mathcal{L}_m)$ theory Eq. (52) and the Einstein-dilaton theory Eq. (53), comes from the fact that any non-linear algebraic functions of the Ricci scalar in the action are equivalent to having an additional scalar degree-of-freedom with gravitational strength [46–48]. To obtain the field equation, we apply the least action principle, ie: $\delta S_{Ed} = 0$. The metric field equation reads

$$ \begin{aligned} R_{\mu\nu} - \frac{1}{\phi}\partial_{\mu}\partial_{\nu}\phi + \frac{1}{\phi}\Gamma^{\lambda}_{\mu\nu}\partial_{\lambda}\phi \\ = \frac{\tilde{\kappa}}{\sqrt{\phi}}\left(T_{\mu\nu} - \frac{1}{2}Tg_{\mu\nu}\right) + \frac{1}{2}\frac{\Box\phi}{\phi}g_{\mu\nu}, \end{aligned} \quad (54) $$

$$ \Box\phi = \frac{\tilde{\kappa}\sqrt{\phi}}{3}(T - \mathcal{L}_m). \quad (55) $$

<footnote n="4">⁴ Given that Mach never formulated a Mach's Principle, there exist many different versions of what Mach's principle could mean [36]. Here we refer to the principle that Einstein used to construct his theory of General Relativity and which he named Mach's principle in [37], see also [38, 39].</footnote>
With $T$ the trace of the stress energy tensor.

Interestingly, even if the action of the theory appears quite different from that of General Relativity, it recovers its phenomenology in many common situations. This comes from what has been called the intrinsic decoupling [49]. This kind of decoupling takes place when one has $\mathcal{L}_m = T$.

Indeed in Eq. (55) if $\mathcal{L}_m = T$, the scalar field is not sourced anymore and the field equation given by Eq. (54) can reduce to those of General Relativity without a cosmological constant. This $\mathcal{L}_m = T$ equality is, for example, a good approximation for a universe made of dust and electromagnetic radiation, just like our universe.

## B. Usual post-Newtonian derivation

The usual post-Newtonian derivation is a pertubative scheme that assumes that the motion of solar system objects are such that $v/c \sim \sqrt{Gm/(rc^2)} \ll 1$ and that $T^{\alpha\beta} = O(c^2, c^1, c^0)$ [50]. It is in particular perfectly relevant in the weak-field limit of the solar-system, but becomes less relevant in the strong-field regime of compact objects such as neutron stars and black-holes. Assuming $\mathcal{L}_m = -\rho$ [17, 51], Eq. (55) becomes

$$ \Box\phi = \tilde{\kappa}\sqrt{\phi}P, \qquad (56) $$

where $P$ is the pressure of the considered fluid that source the gravitational field equations. In either the harmonic or standard post-Newtonian gauges, the first order post-Newtonian metric can be parametrized by a potential $w$, a vector potential $w^i$ and the post-Newtonian parameters $\gamma_{PN}$ and $\beta_{PN}$ [49]

$$ g_{00} = -1 + 2\frac{w}{c^2} - 2\beta_{PN}\frac{w^2}{c^4} + \mathcal{O}(1/c^6), \quad (57a) $$

$$ g_{0i} = -2(1 + \gamma_{PN})\frac{w^i}{c^3} + \mathcal{O}(1/c^5), \quad (57b) $$

$$ g_{ij} = \delta_{ij} \left(1 + 2\gamma_{PN} \frac{w}{c^2}\right) + \mathcal{O}(1/c^4). \quad (57c) $$

Injecting this metric in the field equations (54) and (56) results in

$$ \gamma_{PN} = \beta_{PN} = 1, \quad (58) $$

and

$$ \begin{aligned} w &= w_{GR} - \frac{1}{c^2} G \int \frac{P(\mathbf{x}') d^3 \mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|} + \mathcal{O}(1/c^4), \\ &=: w_{GR} + \frac{1}{c^2} \delta w + \mathcal{O}(1/c^4), \end{aligned} \quad (59a) $$

$$ w^i = w_{GR}^i + \mathcal{O}(1/c^2), \quad (59b) $$

where $w_{GR}$ and $w_{GR}^i$ are the expressions of the potentials predicted by General Relativity. From Eq. (57) and Eq. (59a), one can see that the metric characterizing the considered theory differs from the General Relativity

metric at the first post-Newtonian level by a $1/c^4$ term in the temporal component of the metric. The solution of the scalar-field equation (56) can be written as

$$ \frac{\delta\phi}{\phi_0} = -2\delta w + \mathcal{O}(1/c^2), \quad (60) $$

where $\delta\phi := c^4(\phi - \phi_0)$. In other words, at leading order, one has

$$ A = 1 - \frac{2\omega}{c^2} + O(c^{-4}), \quad (61a) $$

$$ B = 1 + \frac{2\omega}{c^2} + O(c^{-4}), \quad (61b) $$

where

$$ \omega = G \int \frac{\rho(\mathbf{x}') d^3 \mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|} + O(c^{-4}). \quad (62) $$

## C. Deriving an exterior solution of Entangled Relativity stars

The goal of this part is to derive an exterior solution of Entangled Relativity stars using the same methodology as for the Brans-Dicke case explain briefly in Sec II, used by [16]. The logic is exactly the same, only the field equations are slightly changed between the two cases.

Hence we start in Einstein frame with the Janis-Newman-Winicour metric [20] that we re-expressed in the Jordan frame thanks to a conformal transformation. One considers the following rescaling of the scalar field

$$ d\varphi = \frac{\sqrt{3}}{2} d \ln \Phi, \quad (63) $$

while the metric solution reads

$$ ds^2 = -Adt^2 + Bdr^2 + Cd\Omega^2, \quad (64) $$

where the function A, B and C are given by

$$ \left\{ \begin{aligned} A &= -\frac{1}{\Phi} \left(1 - \frac{k}{r}\right)^d, \\ B &= \frac{1}{\Phi} \left(1 - \frac{k}{r}\right)^{-d}, \\ C &= \frac{r^2}{\phi} \left(1 - \frac{k}{r}\right)^{1-d}, \\ \varphi &= -\frac{\sqrt{1-d^2}}{2} \ln \left(1 - \frac{k}{r}\right), \end{aligned} \right. \quad (65) $$

in which $k$ is linked to the mass by $k = \frac{2m}{d}$ and $d$ is the scalar charge. The stationarity and spherical symmetry of the solution imposes an energy-momentum tensor of the form

$$ T_{\nu}^{\mu} = \operatorname{diag}(-\rho, P_{\parallel}, P_{\perp}, P_{\perp}), \quad (66) $$
with $\rho$ the energy density, $P_\parallel$ and $P_\perp$ are respectively the parallel and perpendicular pressure. These are dependent functions of $r$.

To derive the exterior solution, only the (00) component of the field equation Eq.(54) and the scalar field equation Eq.(55) are relevant. They read

$$ \left( \frac{\phi C A_{,r}}{\sqrt{AB}} \right)_{,r} = \frac{\tilde{\kappa} \sqrt{\tilde{\phi}}}{3} (-4T_0^0 + 2T_i^i + \mathcal{L}_m) r^2 \sqrt{AB^3}, \quad (67) $$

$$ \left( \sqrt{\frac{A}{B}} C \phi_{,r} \right)_{,r} = \frac{\tilde{\kappa} \sqrt{\tilde{\phi}}}{3} (T_0^0 + T_i^i - \mathcal{L}_m) r^2 \sqrt{AB^3}, \quad (68) $$

where $X_{,r} := \frac{\partial X}{\partial r}$. As the left-hand side of these equations is in an exact derivative forms, it can easily be integrated from the stars center to a distance $r$ superior to the star's radius $r_*$.

$$ E^* := \int_0^{r^*} dr \frac{\tilde{\kappa} \sqrt{\tilde{\phi}}}{3} (-T_0^0) r^2 \sqrt{AB^3}, \quad (69a) $$

$$ P^* := \int_0^{r^*} dr \frac{\tilde{\kappa} \sqrt{\tilde{\phi}}}{3} T_i^i r^2 \sqrt{AB^3}, \quad (69b) $$

$$ \mathcal{L}_m^* := \int_0^{r^*} dr \frac{\tilde{\kappa} \sqrt{\tilde{\phi}}}{3} \mathcal{L}_m r^2 \sqrt{AB^3}. \quad (69c) $$

Hence Eq.(67,68) becomes

$$ \frac{\phi C A_{,r}}{\sqrt{AB}} = 4E^* + 2P^* + \mathcal{L}_m^*, \quad (70) $$

$$ \sqrt{\frac{A}{B}} C \phi_{,r} = -E^* + P^* - \mathcal{L}_m^*. \quad (71) $$

Now, by solving with the metric Eq.(65), one finds

$$ dk + k \frac{\sqrt{1-d^2}}{\sqrt{3}} = 4E^* + 2P^* + \mathcal{L}_m^*, \quad (72a) $$

$$ k \frac{\sqrt{1-d^2}}{\sqrt{3}} = E^* - P^* + \mathcal{L}_m^*, \quad (72b) $$

and then one finds

$$ k = 3 \sqrt{(E^* + P^*)^2 + \frac{(E^* - P^* + \mathcal{L}_m^*)^2}{3}}, \quad (73a) $$

$$ d = \frac{3(E^* + P^*)}{k}, \quad (73b) $$

as in the Brans-Dicke case one can formulate a parameter $\Theta = \frac{P^*}{E^*}$ and hence rewrite $k$ and $d$ as

$$ k = 3E^* \sqrt{\left(1 + \Theta\right)^2 + \frac{\left(1 - \theta + \frac{\mathcal{L}_m^*}{E^*}\right)^2}{3}}, \quad (74a) $$

$$ d = \frac{1 + \Theta}{\sqrt{\left(1 + \Theta\right)^2 + \frac{\left(1 - \theta + \frac{\mathcal{L}_m^*}{E^*}\right)^2}{3}}}, \quad (74b) $$

Now one can inject Eq. (74) into Eq. (11) as the latter has been obtained through metric expansion in [16].<fnref n="5" /> In the case of Entangled Relativity, $\omega = 0$. Finally one obtains a formula for the exact parameters in Entangled Relativity.

$$ \beta_{\text{exact}} = 1, \quad (75a) $$

$$ \gamma_{\text{exact}} = \frac{1 + 2\Theta - \frac{\mathcal{L}_m^*}{2E^*}}{2 + \Theta + \frac{\mathcal{L}_m^*}{2E^*}}, \quad (75b) $$

$$ \delta_{\text{exact}} = \frac{4}{3} \gamma_{\text{exact}}^2 - \quad (75c) $$

$$ \frac{1}{3} \left( \frac{4 + 2\Theta + \frac{\mathcal{L}_m^*}{E^*}}{\sqrt{9(1 + \Theta)^2 + 3\left(1 - \theta + \frac{\mathcal{L}_m^*}{E^*}\right)^2}} \right)^{-2} . $$

One can especially see that the type of matter will play a crucial role in the resulting curvature because $\mathcal{L}_m^*$ appears in the computation of the exact post-Newtonian parameters. For instance, for an electromagnetic radiation one has $\mathcal{L}_m = T = 0$, whereas for a magnetic field, one has $T \neq \mathcal{L}_m \propto B^2$, where B is the magnetic field. If one considers the example of a dust fluid (i.e. $P = 0$ and $\mathcal{L}_m = -\rho$), one recovers General Relativity with $\gamma_{\text{exact}} = 1$. This is an expected result, since one has $\mathcal{L}_m = -\rho = T$ for a dust fluid, which implies that the scalar field is not sourced in Eq. (55).

## VI. NUMERICAL EVALUATION OF EXACT PARAMETERS IN ENTANGLED RELATIVITY

In this section we will compute the numerical evaluations of ($\gamma_{\text{exact}}$, $\beta_{\text{exact}}$, $\delta_{\text{exact}}$) in the solar system and for neutron stars in Entangled Relativity.

### A. In the solar system

Although the exact post-Newtonian parameters have been derived in the strong-field limit, they can still be

<footnote n="5">⁵ Let us stress once again that the parameters and their analytical expressions are invariant under a change of the radial coordinate—see Sec. II B.</footnote>
derived in the perturbative weak-field case of the usual post-Newtonian developments. In order to derive exact parameters, one needs to know the internal value of A, B and $\phi$, which are given at leading order by Eq. (61) and $\phi = \phi_0 + O(c^{-4})$. Eq. (69a-69b) therefore reduce at leading order to

$$E^* := \frac{\tilde{\kappa}}{3} \int_0^{r_*} \rho r^2 \sqrt{\phi_0} dr + \frac{\tilde{\kappa}}{3} \int_0^{r_*} \rho r^2 \sqrt{\phi_0} \frac{2\omega_{int}}{c^2} dr + O(c^{-4}), \quad (76)$$

$$P^* := \tilde{\kappa} \int_0^{r_*} Pr^2 \sqrt{\phi_0} dr + \tilde{\kappa} \int_0^{r_*} Pr^2 \sqrt{\phi_0} \frac{2\omega_{int}}{c^2} dr + O(c^{-4}), \quad (77)$$

where, assuming a spherical symmetry, one has

$$\omega_{int}(r) = \frac{4\pi G}{r} \int_{0}^{r} \rho(r')r'^{2}dr'. \quad (78)$$

To do the numerical integration, one can use the model S of the sun (⊙) [52] and the PREM model for the earth (⊕) [53], as both gives the pressure and density values of the considered bodies at a considered radius. By computing $\Theta = \frac{P^*}{E^*}$ and using equations (75), one finds, for $\mathcal{L}_m = -\rho$,

$$\gamma_{\text{exact}}^{\odot} - 1 \approx 7.10^{-6}, \quad (79a)$$

$$\gamma_{\text{exact}}^{\oplus} - 1 \approx 9.10^{-10}, \quad (79b)$$

$$\delta_{\text{exact}}^{\odot} - 1 \approx 2.10^{-5}, \quad (79c)$$

$$\delta_{\text{exact}}^{\oplus} - 1 \approx 2.10^{-9}. \quad (79d)$$

The smallness of the values of the exact parameters for Solar System bodies in Entangled Relativity comes from the intrinsic decoupling mentioned in Sec. VA. The codes that we used to compute these values are freely accessible on github at [26]. These predicted values are beyond our current capabilities of measurements, but does not seem to be too far when considering the Sun. Indeed, current constraints are $\gamma_{\odot} = 1+(2.1\pm2.3)\times10^{-5}$ [54]. However, at the $10^{-6}$ level of accuracy for the constraint on $\gamma$, because $\sqrt{Gm_{\odot}/(R_{\odot}c^2)} \sim 10^{-6}$, one can no longer use the leading order development of the propagation of electromagnetic waves and one has to use a full $c^{-4}$ development instead, in order to be fully consistent at the numerical level. This is beyond current best practices [12], although several studies started to address this difficult question already [55-62].

Of course, this is consistent with the perturbative approach of the usual post-Newtonian parameters, as what $\sqrt{Gm_{\odot}/(R_{\odot}c^2)} \sim 10^{-6} \sim \gamma_{\text{exact}}^{\odot} - 1$ means is that $\gamma_{\text{exact}} - 1 = O(c^{-2})$, and therefore can be neglected at the

first post-Newtonian order. This fact changes in the case of ultra-relativistic strong-field compact objects, where $\gamma_{\text{exact}} - 1$ can become somewhat significantly different from zero, as we will know see.

## B. For neutron stars

By Solving the modified Tolman-Oppenheimer-Volkoff equation, one can simulate compact objects in Entangled Relativity. Because it has been shown that near compact object Entangled Relativity predicts possible deviation with respect to General Relativity [17, 19] we decided to model neutron stars. We used a simple polytropic equation of state (as the goal is just to obtain the broad behaviour of $\gamma_{\text{exact}}$ and $\delta_{\text{exact}}$, exactly as explained for the Brans-Dicke case) given by

$$P = K \rho^{\Gamma}, \quad (80)$$

with P the pressure, $\Gamma$ the adiabatic index given by $\Gamma = 5/3$ and a constant $K = 1.475 \times 10^{-3}$ (fm$^3$/MeV)$^{2/3}$ [17, 22]. We conducted the simulation between 100 and 2000 MeV/fm$^3$ but we only kept the simulations that are such that the speed of sound inside the neutron star is less than the conservative limit considered of $c/\sqrt{3}$ [25]. This one corresponds to a central density of 1578 MeV/fm$^3$ for the considered equation of state.

The deviation from unity of $\gamma_{\text{exact}}$ and $\delta_{\text{exact}}$ can be computed with equations (75), and verified as for the Brans-Dicke case with equations (11), with $\omega = 0$. Here again, both parameter has been calculated with the two methods described in sec. IV B, and gives the same results with a maximal relative difference of 0.01% for $\gamma_{\text{exact}}$ and of 0.02% for $\delta_{\text{exact}}$. The resulting plot showing the behavior of $1 - \gamma_{\text{exact}}$ and $1 - \delta_{\text{exact}}$ are presented in Fig. 5 and Fig. 6 in the context of Entangled Relativity.

As we can see in Fig. 5, deviation of $\gamma_{\text{exact}}$ increases with the density. Deviation is observed to range, in absolute value, from 4.7% to 16.1%. Indeed the densest is the body, the higher is the deviation as it can be expected with Eq. (56). The densest is the object the more the scalar field is sourced and hence the more deviation from General Relativity should occur. As expected, concerning the behavior of $\delta_{\text{exact}}$, we can see on Fig. 6 that it follows the same trajectory as for $\gamma_{\text{exact}}$. Deviation is observed to range, in absolute value, from 11.2% to 40.1%.

This results has been obtained without any free parameter but use one assumption: $\mathcal{L}_m = -\rho$. If $\mathcal{L}_m = T$ instead, then no deviation occurs with respect to General Relativity, since in Eq. (55), the scalar field is not sourced and the field equation Eq. (54) becomes the one of General Relativity.

The equation and code to model neutron stars in Entangled Relativity has already been made and computed in [17]. The codes that generates all the figures presented in this manuscript, are freely available on github
[26]. Given that the code integrates the equations of a somewhat recent theory, we use a first integral of the equations in order to further check the integration in Sec. A. We find that the constant from the first integral remains constant at the $10^{-11}$ level, thereby validating the integration at a very good level of accuracy.

<figure><img src="image_6.png" /><figcaption>FIG. 5: Plot of $1 - \gamma_{\text{exact}}$ versus the central density of a Neutron Star. The values of central density ranges from 100 MeV/fm$^3$ to 1578 MeV/fm$^3$. This limit comes from the fact that we want the speed of sound to be less than the conservative limit considered of $c/\sqrt{3}$ [25].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6: Plot of $1 - \delta_{\text{exact}}$ versus the central density of a Neutron Star. The values of central density ranges from 100 MeV/fm$^3$ to 1578 MeV/fm$^3$. This limit comes from the fact that we want the speed of sound to be less than the conservative limit considered of $c/\sqrt{3}$ [25].</figcaption></figure>

## VII. USING DIPOLAR GRAVITATIONAL WAVES TO CONSTRAIN ENTANGLED RELATIVITY

Given the very large values of $1 - \gamma_{\text{exact}}$ predicted by Entangled Relativity for neutron stars when the on-shell matter Lagrangian is assumed to be $\mathcal{L}_m = -\rho$—see Fig. 5—it seems likely that Entangled Relativity could be heavily constrained under this assumption from observations of binary pulsars [63], white dwarf–pulsar systems [33], or triple star systems [21]. Indeed, one would first expect a significant deviation of the Shapiro delay, which is taken into account when fitting the system parameters to pulsar timing data [63]. More importantly, the non-negligible numerical value of $1 - \gamma_{\text{exact}}$ may imply a non-negligible emission of dipolar gravitational radiation in white dwarf–pulsar systems, which is already very tightly constrained [33, 63].

The expectation of significant dipolar emission stems from the fact that such emissions are generically proportional to the square of the Nordtvedt parameter $\eta_{PN} = 4\beta_{PN} - 3\gamma_{PN} - 1$ divided by $1 - \gamma_{PN}$ [33], which vanishes in General Relativity. In Entangled Relativity, one has $\gamma_{PN} = \beta_{PN} = 1$—see Sec. V B—so that $\dot{\eta}_{PN}/(1 - \gamma_{PN}) = 0$ as in General Relativity, and it would therefore seem legitimate not to expect dipolar emission. However, it is likely that dipolar emission from compact objects is sensitive to $\gamma_{\text{exact}}$ and $\beta_{\text{exact}}$, rather than to $\gamma_{PN}$ and $\beta_{PN}$. Hence, one may still expect dipolar emission in white dwarf–pulsar systems in particular.

In this part we will show how to evaluate the variation in orbital period predicted by Entangled Relativity for a binary system composed of a pulsar and a white dwarf emitting dipolar gravitational waves. To accomplish this, we again used our Tolmann-Oppenheimer-Volkoff simulation in Entangled Relativity [17], with the same polytropic equation of state, as well as with another, more realistic, equation of state. Then we used the analytical external solution to deduce the scalar charge as explained in Sec. IV B, as this one intervenes directly in the computation of the variation of the orbital period.

### A. Evaluating the dipolar radiation predicted by Entangled Relativity

With the method in Sec. IV D, one can verify that Eq. (50) also applies to Entangled Relativity, provided that one imposes $\omega_0 = 0$. Following [33], the time variation of the dipolar period $\dot{P}^D$ for a binary system in scalar-tensor theory is given by

$$ \dot{P}^D = -2\pi n_b \mu_* c^{-3} \frac{q}{q+1} \frac{1+e^2/2}{(1-e^2)^{5/2}} (\alpha_P - \alpha_c)^2, \quad (81) $$

where $n_b$ is the orbital frequency, $\mu_*$ is the gravitational parameter, $q$ is the mass ratio between the pulsar and
the companion, *e* is the eccentricity and $\alpha_X \equiv \alpha_{DEF}$ is the Damour-Esposito Farèse scalar charge of the body *X* where *P* stands for *Pulsar* and *c* for *companion*.

We focus on the relativistic Pulsar-White dwarf low eccentricity binary system PSR J1738+0333.<fnref n="6" /> These type of binary system are extremely fruitful for testing alternative theories of gravity. Indeed, they predict a differential in scalar charge that leads to the emission of dipolar gravitational waves. Since the effect of the latter has not been observed yet, these studies can put stringent constraints on scalar-tensor theories of gravity.

<table><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>n<sub>b</sub></td><td>2.05 &times; 10<sup>-4</sup>s<sup>-1</sup></td></tr><tr><td>&mu;<sub>*</sub></td><td>G * M<sub>c</sub> = 2.40 &times; 10<sup>19</sup>m<sup>3</sup>s<sup>-2</sup></td></tr><tr><td>M<sub>p</sub></td><td>computed</td></tr><tr><td>q</td><td>8.1</td></tr><tr><td>e</td><td>3.4 &times; 10<sup>-7</sup></td></tr><tr><td>&alpha;<sub>p</sub></td><td>computed</td></tr><tr><td>&alpha;<sub>c</sub></td><td>0<sup>a</sup></td></tr></tbody></table>

a Using the TOV framework applied to the case of white dwarfs [19] shows that $\alpha_c \in [-2.72 \times 10^{-5}, 1.21 \times 10^{-3}]$. Using these values instead of 0 for $\alpha_c$ leads to negligible corrections to our estimates.

<caption>TABLE I: Sum up of the parameters used in order to compute $\dot{P}^D$ for the system PSR J1738+0333. The parameters come from [33]. Here, $G_*$ is the value of Newton's constant measured thanks to a Cavendish experiment. The pulsar's mass and its scalar charge $\alpha_p$ has been computed in the case of Entangled Relativity while $\alpha_c$ is negligible.</caption>

To numerically evaluate the variation of the orbital period, we use our Tolman-Oppenheimer-Volkoff simulation presented before (see Sec. III or [17]). But because it assumes a too simple polytropic equation of state, we also run the simulation with a more realistic, piece wise, equation of state [64, 65], which has been coded in [66].<fnref n="7" /> Using the values presented in Table I, we ran the simulation for 300 central densities between 100MeV/fm³ and 2000 MeV/fm³ but we only kept simulations for which the speed of sound is less than the conservative limit usually considered of $c/\sqrt{3}$ [25] in the whole neutron star. However, we also limit our analysis to neutron star masses greater than one solar mass.

<figure><img src="image_8.png" /><figcaption>FIG. 7: Plot of $\dot{P}^D$ versus the predicted masses of the Neutron Stars for two different equation of state.</figcaption></figure>

The values obtained from our simulation in Entangled Relativity range from $\dot{P}^D = -1.69 \times 10^{-11} ss^{-1}$ for a mass of $1.40 M_\odot$ to $\dot{P}^D = -1.62 \times 10^{-12} ss^{-1}$ for a mass of $1.01 M_\odot$ for a polytropic equation of state, and $\dot{P}^D = -6.67 \times 10^{-12} ss^{-1}$ for a mass of $1.23 M_\odot$ to $\dot{P}^D = -4.17 \times 10^{-12} ss^{-1}$ for a mass of $1.01 M_\odot$ for a piece wise equation of state. Our low mass estimation comes from the speed-of-sound limit of $c/\sqrt{3}$. Indeed, This limit appears to impose constraints on the masses of the simulated neutron stars. In particular, there is a significant tension—in General Relativity already—between this velocity bound and neutron stars with masses around two solar mass [25]. The details of our simulations are presented in Fig. 7, where we also present the case where one relaxes the speed of sound limit for the piece wise equation of state in order to illustrate how the speed of sound limit impacts the value of the maximal mass attainable.

On the other hand, the data analysis in [33] imply that one has

$$ \dot{P} = (-25.9 \pm 3.2) 10^{-15} ss^{-1}. \qquad (82) $$

Comparing our results presented above and in Fig. 7 to the inferred value Eq. (82) seems to indicate that the case $L_m = -\rho$ should be heavily constrained by this type of systems. But in order to evaluate precisely to what extent, one would have to readjust all the parameters of the binary system within the framework of Entangled Relativity and to also check a wide range of possible equations of state.

<footnote n="6">⁶ In the sense that we use the inferred parameters values of this system in General Relativity to estimate the amount of dipolar gravitational radiation that one could expect for this type of systems. However, one has to keep in mind that all the parameters are model dependent and have to be readjusted within the new theoretical framework.</footnote>

<footnote n="7">⁷ The adapted code to compute $\dot{P}$ from this piece wise equation of state is given in [26].</footnote>
## B. Discussion

At best, our simulations show that two order of magnitudes differ between the predictions of Entangled Relativity with the assumption that $\mathcal{L}_m = -\rho$ (see Fig.7) and the inferred constraint from observational data in Eq. (82). This seems to indicate that the case $\mathcal{L}_m = -\rho$ should be heavily constrained, if not ruled out, by this type of systems. But in order to evaluate precisely to what extent, one would have to readjust all the parameters within the framework of Entangled Relativity and to also check a wide range of possible equations of state. Still, it is unlikely that another set of parameters and equations of state will be able to recover the observed orbital period data—given the orders of magnitude apparent difference.

Even if this result is confirmed, Entangled Relativity would remain viable if $\mathcal{L}_m = T$ instead, as it has sometimes been argued in the literature [67–69]. But then, Entangled Relativity would not predict any deviation from General Relativity for compact objects that are not extremely magnetized. Indeed, in that situation, only the magnetic field would source the scalar-field. Hence a deviation from General Relativity should be expected for magnetars, but the source term of the scalar field is $B^2/2\mu_0$ and is around $10^{11}$ Tesla [70]), which is typically $10^{-7}$ relatively smaller than the typical energy densities of neutron stars—or their pressure, since $\rho \sim 3P/c^2$ for the densest neutron stars. Therefore, a deviation from General Relativity that would entirely be due to a magnetic field with the assumption $\mathcal{L}_m = T$ is expected to be up to seven order of magnitude smaller than one expects from the matter density with the assumption $\mathcal{L}_m = -\rho$, see [19].

## VIII. CONCLUSION

In this work, and for the first time, we have numerically evaluated the so-called *exact* post-Newtonian parameters in scalar–tensor theories of gravity, focusing on Brans–Dicke theory and on Entangled Relativity. Contrary to the usual post-Newtonian parameters, which are perturbative quantities defined in the weak-field limit, the exact parameters are non-perturbative and explicitly depend on the internal structure of the gravitating body through its energy density and pressure. As such, they remain meaningful in the strong-field regime of compact objects, where the standard post-Newtonian expansion becomes insufficient.

We have developed two independent numerical methods to compute these exact parameters: a direct evaluation based on Tolman–Oppenheimer–Volkoff solutions, and an indirect determination through the

matching of the external spacetime to the Janis–Newman–Winicour solution. The excellent agreement between these two approaches provides a strong validation of both the numerical implementation as well as the analytical expressions previously derived in [15, 16].

In the case of Brans–Dicke theory, we have shown that the exact parameters can differ significantly from their usual post-Newtonian counterparts, with relative deviations reaching several tens of percent for sufficiently compact objects. Interestingly, we found that for ultra-relativistic stars the sourcing of the scalar field is reduced due to the relation between pressure and energy density, leading in some cases to smaller deviations from General Relativity than might have been naively expected. These results highlight that strong-field effects cannot be reliably inferred from weak-field post-Newtonian parameters alone.

We then applied the same framework to Entangled Relativity. While the usual post-Newtonian analysis of this theory yields $\gamma_{PN} = \beta_{PN} = 1$, we have shown that this result does not extend to the exact parameters in the strong-field regime. For neutron stars, and under the assumption that the on-shell matter Lagrangian is $\mathcal{L}_m = -\rho$, the deviations of $\gamma_{exact}$ and $\delta_{exact}$ from unity can become substantial, reaching the level of several percent or more, implying indirectly the prediction of dipolar gravitational waves in Entangled Relativity.

In order to probe this different prediction with respect to General Relativity, we showed that the variation in orbital period in Entangled Relativity when assuming $\mathcal{L}_m = -\rho$ ranges from two to three orders of magnitude above the constraint inferred from observational data. This suggests that the theory is likely to be ruled out by existing data, provided that the assumption $\mathcal{L}_m = -\rho$ holds; however, confirming this requires constraining the binary system’s parameters with observational data under various possible equations of state.

If one assumes $\mathcal{L}_m = T$ instead, as advocated in several works, the scalar field is not sourced by ordinary matter and the theory reduces more effectively to General Relativity for compact objects. In that case, observable deviations would only arise for objects endowed with extremely strong electromagnetic fields, such as magnetars, and are in any case expected to be several orders of magnitude smaller than what is found with the assumption $\mathcal{L}_m = -\rho$.

## ACKNOWLEDGEMENT

The authors would like to thank Bertrand Chauvineau for his insightful comments, for suggesting ways to identify the exact moment at which the pressure vanished numerically, and ways for checking the validity of our integration.
[1] C. Brans and R. H. Dicke, “Mach’s principle and a relativistic theory of gravitation,” Phys. Rev. **124**, 925–935 (1961).

[2] Peter G. Bergmann, “Comments on the scalar-tensor theory,” *International Journal of Theoretical Physics* **1**, 25–36 (1968).

[3] Robert V. Wagoner, “Scalar-tensor theory and gravitational waves,” Phys. Rev. D **1**, 3209–3216 (1970).

[4] Kenneth Nordtvedt, “Equivalence principle for massive bodies. ii. theory,” Phys. Rev. **169**, 1017–1025 (1968).

[5] Irwin I. Shapiro, “Fourth test of general relativity,” Phys. Rev. Lett. **13**, 789–791 (1964).

[6] J. H. Taylor and J. M. Weisberg, “A new test of general relativity - Gravitational radiation and the binary pulsar PSR 1913+16,” ApJ **253**, 908–920 (1982).

[7] Abbott (LIGO Scientific Collaboration and Virgo Collaboration), “Observation of gravitational waves from a binary black hole merger,” Phys. Rev. Lett. **116**, 061102 (2016).

[8] The Event Horizon Telescope Collaboration, *The Astrophysical Journal Letters* **875**, L1 (2019).

[9] James G. Williams, Slava G. Turyshev, and Dale H. Boggs, “Lunar Laser Ranging Tests of the Equivalence Principle with the Earth and Moon,” *International Journal of Modern Physics D* **18**, 1129–1175 (2009), arXiv:gr-qc/0507083 [gr-qc].

[10] Clifford M. Will, “The Confrontation between General Relativity and Experiment,” Living Reviews in Relativity **17**, 4 (2014), arXiv:1403.7377 [gr-qc].

[11] LVK collaboration, “Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Reviews in Relativity **23**, 3 (2020).

[12] Agnès Fienga and Olivier Minazzoli, “Testing theories of gravity with planetary ephemerides,” Living Reviews in Relativity **27**, 1 (2024), arXiv:2303.01821 [gr-qc].

[13] Dimitry Ayzenberg, Lindy Blackburn, Richard Brito, Silke Britzen, Avery E. Broderick, Raúl Carballo-Rubio, Vitor Cardoso, Andrew Chael, Koushik Chatterjee, Yi-fan Chen, Pedro V. P. Cunha, Hooman Davoudiasl, Peter B. Denton, Sheperd S. Doeleman, Astrid Eichhorn, Marshall Eubanks, Yun Fang, Arianna Foschi, Christian M. Fromm, Peter Galison, Sushant G. Ghosh, Roman Gold, Leonid I. Gurvits, Shahar Hadar, Aaron Held, Janice Houston, Yichao Hu, Michael D. Johnson, Prashant Kocherlakota, Priyamvada Natarajan, Héctor Olivares, Daniel Palumbo, Dominic W. Pesce, Surjeet Rajendran, Rittick Roy, Saurabh, Lijing Shao, Shammil Tahura, Aditya Tamar, Paul Tiede, Frédéric H. Vincent, Luca Visinelli, Zhiren Wang, Maciek Wielgus, Xiao Xue, Kadri Yakut, Huan Yang, and Ziri Younsi, “Author Correction: Fundamental physics opportunities with the next-generation Event Horizon Telescope,” Living Reviews in Relativity **28**, 7 (2025).

[14] Clifford M. Will, “The confrontation between general relativity and experiment,” Living Reviews in Relativity **17** (2014), 10.12942/lrr-2014-4.

[15] Bertrand Chauvineau and Hoang Ky Nguyen, “The complete exterior spacetime of spherical brans-dicke stars,” Physics Letters B **855**, 138803 (2024).

[16] Bertrand Chauvineau, “Exact post-newtonian parameters about spherical stars in bergmann-wagonernordtvedt scalar-tensor gravity,” Phys. Rev. D **110**, 064060 (2024).

[17] Denis Arruga, Olivier Rousselle, and Olivier Minazzoli, “Compact objects in entangled relativity,” Phys. Rev. D **103**, 024034 (2021), arXiv:2011.14629 [gr-qc].

[18] Denis Arruga and Olivier Minazzoli, “Analytical external spherical solutions in entangled relativity,” *European Physical Journal C* **81**, 1027 (2021), arXiv:2106.03426 [gr-qc].

[19] T. Chehab, O. Minazzoli, and A. Hees, “Variation of planck’s quantum of action in entangled relativity,” *Classical and Quantum Gravity* (2025), 10.1088/1361-6382/ae30c7.

[20] Allen I. Janis, Ezra T. Newman, and Jeffrey Winicour, “Reality of the schwarzschild singularity,” Phys. Rev. Lett. **20**, 878–880 (1968).

[21] G. Voisin, I. Cognard, P. C. C. Freire, N. Wex, L. Guillepot, G. Desvignes, M. Kramer, and G. Theureau, “An improved test of the strong equivalence principle with the pulsar in a triple star system,” *Astronomy and Astrophysics* **638**, A24 (2020), arXiv:2005.01388 [gr-qc].

[22] José D. V. Arbañil, José P. S. Lemos, and Vilson T. Zanchin, “Polytropic spheres with electric charge: Compact stars, the Oppenheimer-Volkoff and Buchdahl limits, and quasiblack holes,” Phys. Rev. D **88**, 084023 (2013), arXiv:1309.4470 [gr-qc].

[23] T Damour and G Esposito-Farese, “Tensor-multi-scalar theories of gravitation,” *Classical and Quantum Gravity* **9**, 2093 (1992).

[24] Olivier Minazzoli, “2PN/RM gauge invariance in Brans-Dicke-like scalar-tensor theories,” *Classical and Quantum Gravity* **29**, 237002 (2012), arXiv:1210.3073 [gr-qc].

[25] Paulo Bedaque and Andrew W. Steiner, “Sound Velocity Bound and Neutron Stars,” Phys. Rev. Lett. **114**, 031103 (2015), arXiv:1408.5116 [nucl-th].

[26] Thomas Chehab and Olivier Minazzoli, “Codes and scripts to compute all the figures,” https://github.com/ThomasChehab/Evaluating_PN (2025).

[27] David Garfinkle, Gary T. Horowitz, and Andrew Strominger, “Charged black holes in string theory,” Phys. Rev. D **43**, 3140–3143 (1991).

[28] Christoph F. E. Holzhey and Frank Wilczek, “Black holes as elementary particles,” *Nuclear Physics B* **380**, 447–477 (1992), arXiv:hep-th/9202014 [hep-th].

[29] Olivier Minazzoli and Edison Santos, “Charged black hole and radiating solutions in entangled relativity,” *European Physical Journal C* **81**, 640 (2021), arXiv:2102.10541 [gr-qc].

[30] Maxime Wavasseur, Théo Abrial, and Olivier Minazzoli, “Slowly rotating and charged black-holes in entangled relativity,” *General Relativity and Gravitation* **57** (2025), 10.1007/s10714-025-03366-5, arXiv:2411.09327 [gr-qc].

[31] Thibault Damour and Gilles Esposito-Farese, “Nonperturbative strong-field effects in tensor-scalar theories of gravitation,” Phys. Rev. Lett. **70**, 2220–2223 (1993).

[32] Thibault Damour and Gilles Esposito-Farèse, “Tensor-scalar gravity and binary-pulsar experiments,” Phys. Rev. D **54**, 1474–1491 (1996), arXiv:gr-qc/9602056 [gr-qc].
[33] Paulo C. C. Freire, Norbert Wex, Gilles Esposito-Farèse, Joris P. W. Verbiest, Matthew Bailes, Bryan A. Jacoby, Michael Kramer, Ingrid H. Stairs, John Antoniadis, and Gemma H. Janssen, “The relativistic pulsar-white dwarf binary PSR J1738+0333 - II. The most stringent test of scalar-tensor gravity,” MNRAS 423, 3328–3343 (2012), arXiv:1205.1450 [astro-ph.GA].

[34] T. Damour and G. Esposito-Farese, “Tensor-multi-scalar theories of gravitation,” Classical and Quantum Gravity 9, 2093–2176 (1992).

[35] Thibault Damour and Gilles Esposito-Farèse, “Tensor-scalar gravity and binary-pulsar experiments,” Physical Review D 54, 1474–1491 (1996).

[36] Mach's Principle: From Newton's Bucket to Quantum Gravity (Birkhäuser, Boston University, 1995).

[37] A. Einstein, “Prinzipielles zur allgemeinen Relativitätstheorie,” Annalen der Physik 360, 241–244 (1918), translation available at https://einsteinpapers.press.princeton.edu/vol7-trans/49.

[38] C. Hoefer, “Einstein’s Formulations of Mach’s Principle,” in Mach’s Principle: From Newton’s Bucket to Quantum Gravity, edited by Julian B. Barbour and Herbert Pfister (Birkhäuser, Boston University, 1995) p. 67.

[39] O. Minazzoli, “On the Principle of Relativity of Inertia in Both General and Entangled Relativities,” Physics of Particles and Nuclei 55, 1488–1493 (2024).

[40] Olivier Minazzoli and Maxime Wavasseur, “Compact objects with scalar charge embedded in a magnetic or electric field in Einstein–Maxwell-dilaton theories,” Eur. Phys. J. C 85, 474 (2025), arXiv:2502.13829 [gr-qc].

[41] Olivier Minazzoli, Maxime Wavasseur, and Thomas Chehab, “Deriving entangled relativity,” Phys. Lett. B 873, 140117 (2026), arXiv:2506.15209 [gr-qc].

[42] Hendrik Ludwig, Olivier Minazzoli, and Salvatore Capozziello, “Merging matter and geometry in the same Lagrangian,” Physics Letters B 751, 576–578 (2015), arXiv:1506.03278 [gr-qc].

[43] Olivier Minazzoli, “Rethinking the link between matter and geometry,” Phys. Rev. D 98, 124020 (2018), arXiv:1811.05845 [gr-qc].

[44] Olivier Minazzoli, “Quantum of action in entangled relativity,” arXiv e-prints , arXiv:2206.03824 (2022), arXiv:2206.03824 [gr-qc].

[45] Olivier Minazzoli, “Rethinking the link between matter and geometry,” Physical Review D 98 (2018), 10.1103/physrevd.98.124020.

[46] S. Capozziello and M. De Laurentis, “F(R) theories of gravitation,” Scholarpedia 10, 31422 (2015), revision #147843.

[47] P. Teyssandier and Ph. Tourrenc, “The Cauchy problem for the $R+R^2$ theories of gravity without torsion,” Journal of Mathematical Physics 24, 2793–2799 (1983).

[48] Andrzej Jakubiec and Jerzy Kijowski, “On theories of gravitation with nonlinear Lagrangians,” Phys. Rev. D 37, 1406–1409 (1988).

[49] Olivier Minazzoli and Aurélien Hees, “Intrinsic Solar System decoupling of a scalar-tensor theory with a universal coupling between the scalar field and the matter Lagrangian,” Phys. Rev. D 88, 041504 (2013), arXiv:1308.2770 [gr-qc].

[50] T. Damour, M. Soffel, and C. Xu, “General-relativistic celestial mechanics. I. Method and definition of reference systems,” Phys. Rev. D 43, 3273–3307 (1991).

[51] Olivier Minazzoli and Tiberiu Harko, “New derivation of the Lagrangian of a perfect fluid with a barotropic equation of state,” Phys. Rev. D 86, 087502 (2012), arXiv:1209.2754 [gr-qc].

[52] J. Christensen-Dalsgaard, W. Däppen, S. V. Ajukov, E. R. Anderson, H. M. Antia, S. Basu, V. A. Baturin, G. Berthomieu, B. Chaboyer, S. M. Chitre, A. N. Cox, P. Demarque, J. Donatowicz, W. A. Dziembowski, M. Gabriel, D. O. Gough, D. B. Guenther, J. A. Guzik, J. W. Harvey, F. Hill, G. Houdek, C. A. Iglesias, A. G. Kosovichev, J. W. Leibacher, P. Morel, C. R. Proffitt, J. Provost, J. Reiter, E. J. Rhodes, F. J. Rogers, I. W. Roxburgh, M. J. Thompson, and R. K. Ulrich, “The current state of solar modeling,” Science 272, 1286–1292 (1996), https://www.science.org/doi/pdf/10.1126/science.272.5266.1286.

[53] Adam M. Dziewonski and Don L. Anderson, “Preliminary reference earth model,” Physics of the Earth and Planetary Interiors 25, 297–356 (1981).

[54] B. Bertotti, L. Iess, and P. Tortora, “A test of general relativity using radio links with the cassini spacecraft,” Nature 425, 374–376 (2003), https://doi.org/10.1038/nature01997.

[55] Neil Ashby and Bruno Bertotti, “Accurate light-time correction due to a gravitating mass,” Classical and Quantum Gravity 27, 145013 (2010), arXiv:0912.2705 [gr-qc].

[56] Olivier Minazzoli and Bertrand Chauvineau, “Scalar-tensor propagation of light in the inner solar system including relevant $c^{-4}$ contributions for ranging and time transfer,” Classical and Quantum Gravity 28, 085010 (2011), arXiv:1007.3942 [gr-qc].

[57] Xue-Mei Deng and Yi Xie, “Two-post-newtonian light propagation in the scalar-tensor theory: An $n$-point mass case,” Phys. Rev. D 86, 044007 (2012).

[58] Bernard Linet and Pierre Teyssandier, “New method for determining the light travel time in static, spherically symmetric spacetimes. Calculation of the terms of order $G^3$,” Classical and Quantum Gravity 30, 175008 (2013), arXiv:1304.3683 [gr-qc].

[59] A. Hees, S. Bertone, and C. Le Poncin-Lafitte, “Relativistic formulation of coordinate light time, Doppler, and astrometric observables up to the second post-Minkowskian order,” Phys. Rev. D 89, 064045 (2014), arXiv:1401.7622 [gr-qc].

[60] B. Linet and P. Teyssandier, “Time transfer functions in Schwarzschild-like metrics in the weak-field limit: A unified description of Shapiro and lensing effects,” Phys. Rev. D 93, 044028 (2016), arXiv:1511.04284 [gr-qc].

[61] Paolo Cappuccio, Ivan di Stefano, Gael Cascioli, and Luciano Iess, “Comparison of light-time formulations in the post-Newtonian framework for the BepiColombo MORE experiment,” Classical and Quantum Gravity 38, 227001 (2021), arXiv:2201.05092 [gr-qc].

[62] Sven Zschocke, “Light propagation in 2PN approximation in the monopole and quadrupole field of a body at rest: Initial value problem,” Phys. Rev. D 105, 024040 (2022), arXiv:2201.06296 [gr-qc].

[63] M. Kramer, I. H. Stairs, R. N. Manchester, N. Wex, A. T. Deller, W. A. Coles, M. Ali, M. Burgay, F. Camilo, I. Cognard, T. Damour, G. Desvignes, R. D. Ferdman, P. C. C. Freire, S. Grondin, L. Guillemot, G. B. Hobbs, G. Janssen, R. Karuppusamy, D. R. Lorimer, A. G. Lyne, J. W. McKee, M. McLaughlin, L. E. Münch, B. B. P.
Perera, N. Pol, A. Possenti, J. Sarkissian, B. W. Stappers, and G. Theureau, “Strong-field gravity tests with the double pulsar,” Phys. Rev. X 11, 041050 (2021).

[64] Potekhin, A. Y., Fantina, A. F., Chamel, N., Pearson, J. M., and Goriely, S., “Analytical representations of unified equations of state for neutron-star matter,” AA 560, A48 (2013).

[65] F. Douchin and P. Haensel, “A unified equation of state of dense matter and neutron star structure,” AA 380, 151–167 (2001), arXiv:astro-ph/0111092 [astro-ph].

[66] Denis Arruga, “Codes and scripts with the piece wise equation of state,” https://github.com/denisArruga/ER_hotspot (2022).

[67] P. P. Avelino and R. P. L. Azevedo, “Perfect fluid Lagrangian and its cosmological implications in theories of gravity with nonminimally coupled matter fields,” Phys. Rev. D 97, 064018 (2018), arXiv:1802.04760 [gr-qc].

[68] P. P. Avelino and R. P. L. Azevedo, “On-shell Lagrangian of an ideal gas,” Phys. Rev. D 105, 104005 (2022), arXiv:2203.04022 [gr-qc].

[69] S. R. Pinto and P. P. Avelino, “Deviations from the von Laue Condition: Implications for the On-Shell Lagrangian of Particles and Fluids,” arXiv e-prints, arXiv:2502.10427 (2025), arXiv:2502.10427 [gr-qc].

[70] Isaac Vidana, “A short walk through the physics of neutron stars,” (2018), arXiv:1805.00837 [nucl-th].

A. VALIDATING OUR INTEGRATION MODEL

Checking the consistency of our integration can be achieved through the use of a first integral. Indeed, by direct integration of Eq. (25), one can obtain a conserved quantity given by

$$C = 5 \ln \left( \frac{c^2}{k^{3/5}} + P^{2/5} \right) + \ln a, \quad (83)$$

where C is a constant. As a reminder, $k = 1.475 \times 10^{-3}$ (fm$^3$/MeV)$^{2/3}$. P is the pressure and $a \equiv g_{tt}$. This quantity can be computed from the initial values and monitored throughout the integration process. Verifying that this constant remains indeed constant allows us to check the consistency of our integration. We find that our integration, both for Brans-Dicke and Entangled Relativity$^8$ theories, leads to relative variations of the constant of order $10^{-11}$ or below$^9$, as illustrated in Fig. 8.

<footnote>$^8$ In Entangled Relativity, Eq. (83) is modified by an additive term $\ln \phi$.</footnote>

<footnote>$^9$ The variation increase with increasing density. A relative variation of the order of $10^{-11}$ is the maximum for both theories.</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 8: Plot of the relative variation of the constant given by $\frac{C_i-C_0}{C_0}$ throughout the numerical integration, versus the radius in Entangled Relativity, for a central density of 100 MeV/fm$^3$.</figcaption></figure>