# A Post-Newtonian Analysis of Regularized 4D-EGB Theory: Complete Set of PPN
Parameters and Observational Constraints

Júnior D. Toniato<EMAIL_ADDRESS>Departamento de Química e Física -
Centro de Ciências Exatas, Naturais e da Saúde, Universidade Federal do
Espírito Santo - Campus Alegre, ES, 29500-000, Brazil.  Núcleo de Astrofísica
e Cosmologia - Cosmo-Ufes, Universidade Federal do Espírito Santo, Vitória,
ES, 29075-910, Brazil.  Martín G. Richarte<EMAIL_ADDRESS>Departamento de
Física - Universidade Federal do Espírito Santo, 29075-910 Vitória, ES, Brazil
PPGCosmo, CCE - Universidade Federal do Espírito Santo, 29075-910 Vitória, ES,
Brazil  Departamento de Física, Facultad de Ciencias Exactas y Naturales,
Universidad de Buenos Aires, Ciudad Universitaria 1428, Pabellón I, Buenos
Aires, Argentina

###### Abstract

We performed a post-Newtonian analysis of the regularized four-dimensional
Einstein-Gauss-Bonnet gravitational theory (4D-EGB). The resulting metric
differs from the classical parametrized post-Newtonian (PPN) formalism in that
a new gravitational potential arises from the integration of the approximate
field equations. We also investigated the conserved quantities and equations
of motion for massive bodies and light rays to a certain degree. By computing
the predicted periastron advance rate in a binary system, we obtained an
observational constraint that is stronger than those of previous analyses.
Although the usual 10 PPN parameters can still be derived within the PPN
framework, an extra parameter is needed to account for the full post-Newtonian
tests.

##  I Introduction

Einstein’s theory of general relativity (GR) has passed numerous experimental
tests since its inception in the early 20th century. GR could accurately
explain the precession of Mercury’s perihelia, a significant accomplishment at
the time. Over the past century, GR has been put to the test on various fronts
[1, 2]. In 2015, a major milestone was achieved when the LIGO interferometer
successfully detected the first signals of gravitational waves (GWs) emitted
by the merger of two black holes [3]. This groundbreaking discovery was
followed by another significant event in 2017, when Advanced LIGO and Virgo
detected the electromagnetic counterpart of a binary neutron star coalescence
[4]. In 2019, the EHT collaboration released the first image of the shadow
cast by the supermassive black hole harbored by the elliptical galaxy M87
belonging to the nearby Virgo galaxy cluster [5]. These remarkable
achievements have opened up new avenues for studying the universe and have
further validated the predictions of Einstein’s theory of general relativity.

However, there are several areas where (GR) still requires improvement in
order to accurately describe various physical phenomena. One such area
pertains to the observed acceleration of the universe, as indicated by the
seminal works of type Ia supernovae (standard candles) [6, 7]. This suggests
the existence of a new component (dark energy) with negative pressure, which
is responsible for this acceleration. Initially, the addition of a
cosmological constant to Einstein’s field equation was seen as the most
plausible solution, despite the discrepancy between the observed value and the
theoretical prediction from vacuum energy in QFT [8, 9, 10, 11]. Another
avenues consider that dark energy could be explained in terms of a dynamical
field with really a small mass [12, 13, 14]. A modified theory of gravity
could also explain the universe’s acceleration on a cosmological scale (see
[15, 16, 17], and [18]). A more reliable approach could be introducing a new
degree of freedom, a scalar field, coupled to curvature terms without spoiling
all the nice properties of Einstein’s field equations (see the conditions
mentioned in Lovelock’s theorem [19]) and without introducing Ostrogradsky
instability [20, 21]. Among the new healthy theories with all those
properties, the so-called Hordensky’s theories [22] have attracted
considerable attention in recent years [23, 24, 25, 26]. From an observational
standpoint, there are compelling reasons to consider theories such as
Hordensky’s and beyond. These theories have successfully met rigorous
constraints based on multi-messenger gravitational-wave astronomy. The
detection of the GW170817 binary neutron star merger, along with its
associated electromagnetic counterparts, has provided an accurate bound on the
speed of gravitational waves, $|c_{g}/c-1|<10^{-16}$ [27, 28].

In connection with the frameworks mentioned above, an appealing version of the
Einstein-Gauss-Bonnet theory in four dimensions (4D-EGB) was constructed after
a regularization process [29, 30, 31, 32]. As is well-known, the GB curvature
invariant reduces to a topological surface term in four dimensions, so whether
the EGB theory can be defined correctly in the limit $D\rightarrow 4$ was an
open question [33, 34, 35, 36, 37]. A neat way to deal with this issue is by
considering a conformal transformation and employing a subtraction method to
render the theory well-defined in that limit [37]. The conformal factor
involves a scalar field that introduces a new degree of freedom into the
model. Interestingly enough, the regularized 4D-EGB model can be accommodated
as a particular case within the shift-symmetric Horndeski’s theories [22, 26].
Another possibility for obtaining a well-defined 4D-EGB theory can be achieved
by implementing a Kaluza-Klein dimensional reduction from a higher-dimensional
Einstein-Gauss-Bonnet gravity plus a scalar field. In that case, the scalar
field encodes the information of the volume’s size of the compact internal
dimension with maximal symmetry [31].

One way to characterize new theories is by using observational tests on
different scales, from the solar system to the cosmological scales [38]. In
that regard, the regularized 4D-EGB model was scrutinized on several fronts.
Considering the limits on the GW velocity from the observations of GW170817
and GRB 170817A, the bounds on the rescaled coupling parameter reads
$0<\hat{\alpha}<10^{50}\,{\rm{eV}}^{-2}$ [39]. Stronger limits on positive
values come from binary black hole systems, resulting in the overall range
$\hat{\alpha}<10^{8}\,{\rm{m}}^{2}$ [38]. In the weak-field regime where the
gravitational field is considerably tiny, using the LAGEOS satellite data, one
can obtain comparable bounds on the Gauss-Bonnet coupling,
$\hat{\alpha}<10^{10}{\rm{m}}^{2}$ [38]. Besides, assuming that the lightest
component of the GW190814 event corresponds to black hole with mass
$M=2.59^{+0.08}_{-0.09}\,M_{\odot}$, then the bound on GB coupling reduces
several orders of magnitude, $\hat{\alpha}<59\times 10^{6}{\rm{m}}^{2}$
[40].111The most stringent limit on negative value of GB coupling comes from
the condition that atomic nuclei must not be shielded by a horizon, yielding
$-\hat{\alpha}>10^{-83}H^{2}_{0}$ [40]. Besides, primordial GWs also provide
another similar constraint on the GB coupling,
$|\hat{\alpha}|<10^{36}{\rm{m}}^{2}$ [41],[37]

In 2024, new constraints were reported within the context of FRW background
with non-zero spatial curvature [42]. To be more precise, the latter case
leads to a modified Friedmann equation for the Hubble parameter after
integrating the scalar field dynamical equation; in fact, the non-linear
Friedmann equation has a new dark radiation term multiplied by the GB coupling
constant. Physically speaking, the extra-dark radiation term will modify the
(massless) species number of neutrinos at both the background and perturbative
levels [42]. The constraints based on the ACTPol data alone implies
$\hat{\alpha}C^{2}=(-9\pm 6)\times 10^{-6}H^{2}_{0}$, where
$\hat{\alpha}C^{2}=\Omega^{\rm{extra}}_{\rm{rad}}H^{2}_{0}$ stands for the
geometrical dark radiation term [42].

We aim to extend previous works on astrophysical constraints on the
regularized four-dimensional Einstein-Gauss-Bonnet theory by reevaluating the
parametrized post-Newtonian (PPN) formalism [43] within the context of the
4D-EGB theory and adding stronger astrophysical observational constraints
along with a whole discussion about bounds comings from complementary probes.

The paper proceeds as follows. We obtain the effective field equations for the
metric and scalar fields in Section II. Section III is devoted to the
determination of the post-Newtonian form of the metric within the PPN
formalism standards, where conclusions on the $\gamma$ parameter are already
obtained. We explore in Section IV the post-Newtonian conservation laws
statements to drawn conclusions on the $\zeta$’s and $\alpha$’s parameters.
All those parameters does maintain their physical meaning unaltered and it is
shown to have the same value as in GR. Section V is dedicated to computing the
center-of-mass acceleration of an extended object in a system of $N$-bodies
within the PN approximation, where we capture the standard GR term plus
additional terms proportional to the Gauss-Bonnet coupling, $\hat{\alpha}$.
With the previous results, we can estimate the shift in the orbital periastron
for a binary system due to the GB coupling, which is done in Section VI. We
also examine the periastron advance rate in the cases of Mercury and the
double pulsar. Section VII discusses the role played by the $\beta$ parameter
in the PPN approach and the Nordtvedt effect, along with an analysis on haw to
extend th formalism to include the corrections of the EGB coupling. We state
our conclusions in Section VIII.

##  II The regularized 4D-EGB theory

In this section, we will derive the field equations of the regularized 4D-EGB
theory as obtained in the previous Refs.[29, 30, 31, 32]. Our starting point
is to consider the gravity sector of the EGB without cosmological constant in
$D$ dimensions,

$S=\int_{\mathcal{M}}d^{D}x\frac{\sqrt{-g}}{2\kappa}\left(R+\alpha\mathcal{G}\right)\
,$ (1)

where $\kappa=8\pi G/c^{4}$, with $G$ the gravitational constant and $c$ the
velocity of light in vacuum whereas $\alpha$ stands for a coupling constant
associated to the quadratic term, $\mathcal{G}$. The regularization process
involves a conformal transformation, followed by a redefinition of the
coupling constant. Specifically, we apply
$\tilde{g}_{\mu\nu}=e^{2\Phi}g_{\mu\nu}$, which then transforms the square
root of the determinant as $\sqrt{-\tilde{g}}=e^{D\Phi}\sqrt{-g}$. The Gauss-
Bonnet term changes as follows [37],

$\displaystyle\sqrt{-\tilde{g}}\tilde{\mathcal{G}}$
$\displaystyle=e^{(D-4)\Phi}\sqrt{-g}\Big{[}\mathcal{G}-8(D-3)R^{\mu\nu}(\nabla_{\mu}\Phi\nabla_{\nu}\Phi-\nabla_{\mu}\nabla_{\nu}\Phi)-2(D-3)(D-4)R(\nabla\Phi)^{2}$
$\displaystyle+4(D-2)(D-3)^{2}\Box\Phi(\nabla\Phi)^{2}-4(D-2)(D-3)\nabla_{\mu}\nabla_{\nu}\Phi\nabla^{\mu}\nabla^{\nu}\Phi-4(D-3)R\Box\Phi$
$\displaystyle+4(D-2)(D-3)(\Box\Phi)^{2}+8(D-3)(\nabla_{\mu}\Phi\nabla_{\nu}\Phi)(\nabla_{\mu}\nabla_{\nu}\Phi)+(D-1)(D-2)(D-3)(D-4)(\nabla\Phi)^{4}\Big{]},$
(2)

The regularization method requires to subtract a counter-term
$\alpha\sqrt{-\tilde{g}}\tilde{\mathcal{G}}$ in the original action. The idea
is to expand the exponential around $D=4$ and keep the lowest order in $(D-4)$
and neglect all higher order terms, namely,
$e^{(D-4)\Phi}=1+(D-4)\Phi+\mathcal{O}[(D-4)^{2}]$. We use the Bianchi
identities along with the Bochner’s formula in curved spacetime,
$\frac{1}{2}\nabla^{\mu}\left(\nabla\Phi\right)^{2}=\left(\nabla_{\mu}\nabla_{\nu}\Phi\right)^{2}+\nabla_{\mu}\Big{(}\Box\Phi\nabla^{\mu}\Phi\Big{)}-\left(\Box\Phi\right)^{2}+R^{\mu\nu}\nabla_{\mu}\Phi\nabla_{\nu}\Phi$,
for simplifying the action. In doing so, we integrate by parts, discarding any
surface terms. We set $D\rightarrow 4$ for the terms that are not of the form
$D-4$, and the divergences in the original action is canceled provided the
term proportional to $(D-4)$ remains finite as they are multiplied by the re-
scaled coupling $\hat{\alpha}/(D-4)$ [29, 30, 31, 32]. The resulting action of
the system reads [29, 30, 31, 32, 44],

$\displaystyle S^{\text{4D-EGB}}_{\text{reg}}=\int
d^{4}x\frac{\sqrt{-g}}{2\kappa}\big{[}$ $\displaystyle
R+\hat{\alpha}\big{(}4G^{\mu\nu}\nabla_{\mu}\Phi\nabla_{\nu}\Phi-\Phi\mathcal{G}$
$\displaystyle+4\square\Phi(\nabla\Phi)^{2}+2(\nabla\Phi)^{4}\big{)}\big{]}+S_{m},$
(3)

where we identify $G_{\mu\nu}$ as the Einstein tensor whereas ${\cal
G}=R_{\alpha\beta\mu\nu}R^{\alpha\beta\mu\nu}-4R_{\mu\nu}R^{\mu\nu}+R^{2}$
denotes the Gauss-Bonnet term. As usual, $R$ refers to the Ricci scalar,
$S_{m}$ describes the matter action, and $\nabla_{\mu}$ indicates a covariant
derivative. The extra degree of freedom,$\Phi$, comes from the scalar field
that encodes the conformal transformations that allow us to regularize the
gravitational action for $D=4$ dimensions [29]. For an entire debate on the
nature of this regularization process, we refer the reader to the following
works [29, 30, 31, 32, 44]. To explore the new modifications introduced by
this regularization process, we need to obtain the field’s equation of the
regularized 4D-EGB theory (3). We carry on by computing the variation of
action (3) with respect to the metric, it is given by

$G_{\mu\nu}=\hat{\alpha}\hat{\mathcal{H}}_{\mu\nu}+\kappa T_{\mu\nu}\,,$ (4)

where $T_{\mu\nu}$ is the energy-momentum tensor associated with the matter
fields and the tensor $\hat{\mathcal{H}}_{\mu\nu}$ can be written as

$\displaystyle\hat{\mathcal{H}}_{\mu\nu}=$ $\displaystyle\
2R(\nabla_{\mu}\nabla_{\nu}\Phi-\nabla_{\mu}\Phi\nabla_{\nu}\Phi)+2G_{\mu\nu}\left(\left(\nabla\Phi\right)^{2}-2\Box\Phi\right)+4G_{\nu\alpha}\left(\nabla^{\alpha}\nabla_{\mu}\Phi-\nabla^{\alpha}\Phi\nabla_{\mu}\Phi\right)$
(5)
$\displaystyle+4G_{\mu\alpha}\left(\nabla^{\alpha}\nabla_{\nu}\Phi-\nabla^{\alpha}\Phi\nabla_{\nu}\Phi\right)+4R_{\mu\alpha\nu\beta}\left(\nabla^{\beta}\nabla^{\alpha}\Phi-\nabla^{\alpha}\Phi\nabla^{\beta}\Phi\right)+4\nabla_{\alpha}\nabla_{\nu}\Phi\left(\nabla^{\alpha}\Phi\nabla_{\mu}\Phi-\nabla^{\alpha}\nabla_{\mu}\Phi\right)$
$\displaystyle+4\nabla_{\alpha}\nabla_{\mu}\Phi\nabla^{\alpha}\Phi\nabla_{\nu}\Phi-4\nabla_{\mu}\Phi\nabla_{\nu}\Phi\left(\left(\nabla\Phi\right)^{2}+\Box\Phi\right)+4\Box\Phi\nabla_{\nu}\nabla_{\mu}\Phi-
g_{\mu\nu}\Big{(}2R\left(\Box\Phi-\left(\nabla\Phi\right)^{2}\right)$
$\displaystyle+4G^{\alpha\beta}\left(\nabla_{\beta}\nabla_{\alpha}\Phi-\nabla_{\alpha}\Phi\nabla_{\beta}\Phi\right)+2(\Box\Phi)^{2}-\left(\nabla\Phi\right)^{4}+2\nabla_{\beta}\nabla_{\alpha}\Phi\left(2\nabla^{\alpha}\Phi\nabla^{\beta}\Phi-\nabla^{\beta}\nabla^{\alpha}\Phi\right)\Big{)}.$

The variation of the action (3) with respect to the scalar field, $\Phi$,
leads to its equation of motion,

$\displaystyle
R^{\mu\nu}\nabla_{\mu}\Phi\nabla_{\nu}\Phi-G^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\Phi-\Box\Phi\left(\nabla\Phi\right)^{2}+\nabla_{\mu}\nabla_{\nu}\Phi\nabla^{\mu}\nabla^{\nu}\Phi-(\Box\Phi)^{2}-2\nabla_{\mu}\Phi\nabla_{\nu}\Phi\nabla^{\mu}\nabla^{\nu}\Phi=\frac{1}{8}\mathcal{G}.$
(6)

The field equation (4) can be rearranged with the help of its trace by writing
the Ricci scalar in terms of $T$ and $\hat{\mathcal{H}}$, i.e., the respective
traces of the energy-momentum tensor and the tensor
$\hat{\mathcal{H}}_{\mu\nu}$. This leads us to recast the metric field
equations in a more suitable way for solving them within the post-Newtonian
approximation as the standard Einstein’s field equation plus a new
contribution encoded in the $\hat{\mathcal{H}}_{\mu\nu}$ tensor,

$R_{\mu\nu}=8\pi\left(T_{\mu\nu}-\frac{T}{2}\,g_{\mu\nu}\right)+\hat{\alpha}\left(\hat{\mathcal{H}}_{\mu\nu}-\frac{\hat{\mathcal{H}}}{2}\,g_{\mu\nu}\right),$
(7)

where we chose a geometrodynamics system of units by selecting $G=c=1$. The
aforesaid theory (7) is a particular case of a broader class of shift-
symmetric Horndeski’s theories. These theories have the property that their
equations of motion do not change when a constant shifts the scalar field,
$\Phi\rightarrow\Phi+c$, where $c$ is a constant; therefore, the theory admits
a conserved Noether current [22, 26]. Specifically, it corresponds to
$G_{2}=8\hat{\alpha}X^{2}$, $G_{3}=8\hat{\alpha}X$, $G_{4}=1+4\hat{\alpha}X$,
and $G_{5}=4\hat{\alpha}\ln X$ (where
$X=-\frac{1}{2}\nabla_{\mu}\Phi\nabla^{\mu}\Phi$). General Horndeski theories
have been previously analyzed within the post-Newtonian approximation scheme,
but only by considering models in which the $G_{n}$ functions can be expressed
as a Taylor series around $X=0$; this is not the case for $G_{5}$ [45].
Notwithstanding, equations (4) and (6) can be regularly treated within the
post-Newtonian limit, as seen in the following sections.

##  III Post-Newtonian expansion

Let us focus on solving the theory’s field equation in a regime where the
gravitational field is weak everywhere and matter moves in small velocities
compared to light speed. We proceed by expanding equations (7) and (6) using
the following perturbation’s scheme,

$\displaystyle g_{\mu\nu}$ $\displaystyle=\eta_{\mu\nu}+h_{\mu\nu}\,,$ (8)
$\displaystyle\Phi$ $\displaystyle=\phi_{0}+\phi,$ (9)

where $h_{\mu\nu}$ encodes the metric perturbations concerning Minkowski’s
background, whereas $\phi$ stands for the scalar field perturbation around its
constant background value, $\phi_{0}$.

For the matter content, we use a perfect fluid approach,

$\displaystyle
T^{\mu\nu}=\left(\rho+\rho\Pi+p\right)u^{\mu}u^{\nu}+pg^{\mu\nu}\,,$ (10)

where $\rho$ is the mass density, $\Pi$ is the fluid’s internal energy per
unity mass, $p$ is the pressure and $u^{\mu}=u^{0}(1,\bm{v})$ is the four-
velocity of the fluid. In the post-Newtonian approximation scheme, the energy-
momentum tensor components should be expanded in power of $v$, taking into
account also that $p/\rho\sim\Pi\sim O(2)$. Here and after, we use the
notation $O(N)$ to represent quantities of order $v^{N}$.Moreover, time
derivatives must be considered as order of magnitude one, $\partial_{t}\sim
O(1)$, since the dynamical time scale in the solar system is governed by the
motion of planets.

As expected, the metric perturbations $h_{\mu\nu}$ will also be expanded in
powers of $v$. To obtain the first post-Newtonian corrections to the equations
of motion of massive bodies, one needs to know $h_{00}$ up to order $v^{4}$,
$h_{0i}$ up to order $v^{3}$ and $h_{ij}$ up to order $v^{2}$, where Latin
indices run over the three spatial dimensions. The approximation procedure
described above is standard in the post-Newtonian approach and is explained in
great detail in Refs. [43, 46], which we will follow closely.

The fluid dynamics is subjected not only to conservation of the energy-
momentum tensor, $\nabla_{\nu}T^{\mu\nu}=0$, but also to the conservation of
rest mass density,

$\nabla_{\mu}(\rho u^{\mu})=0\,.$ (11)

This equation can be re-expressed as an effective flat-space continuity
equation as follows,

$\partial_{t}(\rho^{*})+\partial_{i}(\rho^{*}v^{i})=0,$ (12)

where the conserved density is defined in terms of zero component of the four-
velocity field and the square root of the determinant,

$\rho^{*}\equiv{u^{0}}\sqrt{-g}\,\rho.$ (13)

Since $\rho^{*}$ is the conserved density, in the sense of (12), it is more
convenient to use it to express the energy-momentum tensor components. As will
become clear later, having the gravitational potentials expressed in terms of
the conserved density will lead to an easy way to integrate the equations of
motion of massive bodies.

Applying the above approximations, one can verify that the first non-vanishing
terms in equation (6) and for $\hat{\mathcal{H}}_{\mu\nu}$ occur in the fourth
order only. Consequently, 4D-EGB gravity will deviates from GR only at the
fourth order [see Eq. (7)]. Thus, up to $O(3)$, the metric components reads,

$\displaystyle g_{00}$ $\displaystyle=-1+2U+O(4),$ (14) $\displaystyle g_{0i}$
$\displaystyle=-4\,V^{i}+O(5),$ (15) $\displaystyle g_{ij}$
$\displaystyle=(1+2U)\delta_{ij}+O(4),$ (16)

where $U$ is the negative of the Newtonian potential,

$U=\int\frac{\rho^{*}(t,\bm{x}^{\prime})\,d^{3}x^{\prime}}{|\bm{x}-\bm{x}^{\prime}|},$
(17)

and,

$V^{i}=\int\frac{\rho^{*}(t,\bm{x}^{\prime})v^{i}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime}.$
(18)

To determine the fourth-order term of $g_{00}$ we first consider equation (6),
which gives,

$\partial_{ij}\phi\,\partial^{ij}\phi-\left(\nabla^{2}\phi\right)^{2}=\partial_{ij}U\,\partial^{ij}U-\left(\nabla^{2}U\right)^{2},$
(19)

where we are using the notation $\partial_{ij}=\partial/\partial x^{i}\partial
x^{j}$. The simplest solution to the above equations reads,

$\phi=\pm\,U.$ (20)

Using Eqs. (14)-(16) and (20), the $00$-component of field equation (7)
becomes,

$\displaystyle\nabla^{2}h_{00}=-\rho^{*}\left(1+\frac{3}{2}v^{2}-U+\Pi+3p\right)$
$\displaystyle-2\partial_{tt}U\,\pm\,\hat{\alpha}\,\mathcal{G}^{(4)}+O(6),$
(21)

where $\mathcal{G}^{(4)}$ is the fourth-order approximation of Gauss-Bonnet
invariant,

$\mathcal{G}^{(4)}=8\,\partial_{ij}U\,\partial^{ij}U-8\left(\nabla^{2}U\right)^{2}.$
(22)

Assuming asymptotically flatness, the general solution of (21) reads,

$h_{00}=-2U+2(\Psi-U^{2}\mp\hat{\alpha}\,\Phi_{\mathcal{G}})+O(6),$ (23)

with

$\Psi\equiv\dfrac{3}{2}\,\Phi_{1}-\Phi_{2}+\Phi_{3}+3\Phi_{4}+\frac{1}{2}\partial_{tt}X,\\\
$ (24)

and the PN potentials that appear can be recast as follows,

$\displaystyle\Phi_{1}=\int\dfrac{\rho^{*}{}^{\prime}v^{\prime
2}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime},\quad\Phi_{2}=\int\dfrac{\rho^{*}{}^{\prime}U^{\prime}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime},$
(25)
$\displaystyle\Phi_{3}=\int\dfrac{\rho^{*}{}^{\prime}\Pi^{\prime}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime},\quad\Phi_{4}=\int\dfrac{p^{\prime}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime}\,,$
(26) $\displaystyle
X=\int\rho^{*}{}^{\prime}|\bm{x}-\bm{x}^{\prime}|\,d^{3}x^{\prime},\quad\Phi_{\mathcal{G}}=\frac{1}{4\pi}\int\dfrac{\mathcal{G}^{(4)}\,d^{3}x^{\prime}}{|\bm{x}-\bm{x}^{\prime}|}\,.$
(27)

The prime symbol indicates that a quantity is evaluated at the spatial
coordinates $x^{\prime}$.

In brief, the post-Newtonian metric of 4D-EGB theory is given by

$\displaystyle g_{00}$
$\displaystyle=-1+2U+2(\Psi-U^{2})\mp\hat{\alpha}\,\Phi_{\mathcal{G}}+O(6),$
(28) $\displaystyle g_{0i}$ $\displaystyle=-4\,V^{i}+O(5),$ (29)
$\displaystyle g_{ij}$ $\displaystyle=(1+2U)\delta_{ij}+O(4).$ (30)

The above metric is consistent with previous analyses on weak-field and slow-
motion approximations of 4D-EGB gravity [38]. However, several differences
deserve to be mentioned: $i)$ we are working with harmonic gauge instead of
the traditional PPN gauge, and $ii)$ we are defining potentials in terms of
the conserved density. Notwithstanding, this does not change the form of the
correction brought by this theory, represented by the potential
$\Phi_{\mathcal{G}}$. This potential is not covered by the PPN formalism,
implying that it is not correct to read off the PPN parameters directly from
the metric structure. One must investigate the equations of motion of massive
bodies and light rays to correctly identify the new potential influence in the
standard PPN parameters.

It is worth to mention this is not a particularity of 4D-EGB model, once
several modern theories of gravity does not fit entirely to the PPN formalism
assumptions. It seems not to exist a fundamental rule that explains or
determines when a model may include extra-PPN potentials on it. An a priori
simple theory, like purely scalar gravity, can produces an extremely knotty PN
metric [47], while GR extensions motivated by nontrivial renormalization group
effects at large scales may not spoil the PPN framework at all [48].
Determinations of PPN parameters for theories with extra-PPN potentials is
usually possible, as in the case of Palatini $f(R)$ gravity [49], for
instance. But sometimes some extra although reasonable assumptions are needed,
as in scalar-tensor and metric $f(R)$ formulations [50, 51]. Fortunately, the
latter is not the case of 4D-EGB.

###  III.1 The $\bm{\gamma}$ parameter

For the propagation of light, the first post-Newtonian correction is obtained
by considering only the second-order terms of the metric. The PPN formalism
only assumes the Newtonian potential to be present at order two, and any
deviation from GR is encoded in the single parameter $\gamma$ multiplying the
$g_{ii}$ components (see Ref. [50] for a detailed discussion on the
parametrization of propagation of light). Thus, once the 4D-EGB theory
produces only the Newtonian potential at the second-order metric, one can
conclude that the physical meaning of the $\gamma$ parameter remains
unaltered. Since the spatial component of the second-order metric is the same
as in GR [eq. (30)], it is obtained

$\gamma=1,$ (31)

and a full agreement with tests in the solar system based on propagation of
light. The remaining PPN parameters can be inferred from the calculations of
the conserved quantities and the equations of motion of massive bodies. In the
following sections, we will perform this analysis.

##  IV Conserved Quantities

In this section, we obtain the expressions for a fluid system’s total mass-
energy and momentum, which are conserved in the post-Newtonian limit. The re-
scaled density $\rho^{*}$ allows one to define the fluid material mass in a
given volume as follows,

$m=\int\rho^{*}d^{3}x\,.$ (32)

Hence, for an isolated body, by using the continuity equation for $\rho^{*}$
[cf. eq. (12)], it is shown that mass $m$ is conserved, i.e.

$\dfrac{dm}{dt}=0\,.$ (33)

The above result assumes that density $\rho^{*}$ is zero at the boundary of
the integration volume. The statement above only depends on the fluid
definition, not on the used gravity model. The situation is different for
energy and momentum conservation.

From the energy-momentum tensor conservation, we write

$\nabla_{\nu}T^{\mu\nu}=\partial_{\nu}(\sqrt{-g}\,T^{\mu\nu})+\Gamma^{\mu}_{\lambda\nu}(\sqrt{-g}\,T^{\lambda\nu})=0.$
(34)

We want to integrate the above expression in its post-Newtonian approximation
and obtain the conserved quantities in 4D-EGB gravity. For practical reasons,
one can work with a redefinition of eq. (28) such as,

$g_{00}=(g_{00})_{\scriptscriptstyle\rm
GR}\mp\hat{\alpha}\,\Phi_{\mathcal{G}}/2,$ (35)

where $(g_{00})_{\scriptscriptstyle\rm GR}$ stands for the post-Newtonian
metric from GR. Then, one can easily separate the contributions inserted by
the 4D-EGB theory and use the standard GR results as shown in Ref. [46] for
instance.

The statement of energy conservation is derived from the $\mu=0$ component of
(34), where there is no fourth-order metric contributions. Hence the results
are the same as in GR: the leading order reproduces the continuity equation
(12) and the next order gives the energy which is conserved at post-Newtonian
regime, namely

$E\equiv\int\left(\frac{1}{2}\rho^{*}v^{2}+\rho^{*}\Pi-\frac{1}{2}\rho^{*}U\right)d^{3}x.$
(36)

Thus, the total mass-energy is given by

$M=m+E.$ (37)

For the conservation of the total momentum, we consider only the case $\mu=i$
in eq. (34), which can be written as follows,

$(\nabla_{\nu}T^{\nu}_{~{}i})_{\scriptscriptstyle\rm
GR}\pm\dfrac{\hat{\alpha}}{2}\,\rho^{*}\nabla_{i}\Phi_{\mathcal{G}}=0,$ (38)

where $(\nabla_{\nu}T^{\nu}_{~{}i})_{\scriptscriptstyle\rm GR}$ stands, once
again, for the ordinary GR terms. By integrating this expression over the
volume occupied by the fluid, we are able to derive the total momentum which
is conserved. We only need to integrate the last term in (38) since the
results from the GR contributions are already known [46]. To do so, we first
verify that (see Ref.[38]) $\Phi_{\mathcal{G}}$ can be expressed in terms of
$U$ and $\Phi_{7}$ as follows,

$\Phi_{\mathcal{G}}=-8\left(\frac{1}{2}|\nabla U|^{2}-\Phi_{7}\right)\,,$ (39)

where we have defined,222Symbols $\Phi_{5}$ and $\Phi_{6}$ are already used
within GR to represent other gravitational potentials. This is the reason why
we use $\Phi_{7}$ to represent this new potential. It is the same potential
described in Ref. [38] as $\psi_{1}$.

$\Phi_{7}=-\int\rho^{*}{}^{\prime}\partial^{\prime}_{j}U^{\prime}\frac{({x}-{x}^{\prime})^{j}}{|{\bf
x}-{\bf x}^{\prime}|^{3}}d^{3}x^{\prime}.$ (40)

By integrating the gradient of $\Phi_{7}$ and using the identity
$\partial_{i}f(x-x^{\prime})=-\partial_{i}^{\prime}f(x-x^{\prime})$ plus an
exchange of variables, $x\leftrightarrow x^{\prime}$, one arrives at the
following expression,

$\displaystyle\int\rho^{*}\partial_{i}\Phi_{7}\,d^{3}x$
$\displaystyle=\int\rho^{*}\partial^{j}U\partial_{ij}U\,d^{3}x,$
$\displaystyle=\frac{1}{2}\int\rho^{*}\partial_{i}(|\nabla U|^{2})d^{3}x.$
(41)

Combining (39) and (41), one shows that the following relation holds true,

$\int\rho^{*}\partial_{i}\Phi_{\mathcal{G}}\,d^{3}x=0.$ (42)

Consequently, after integrating (38), one gets the same results as in GR for
the total conserved momentum,

$\displaystyle P^{j}=$
$\displaystyle\int\rho^{*}v^{j}\bigg{(}1+\frac{v^{2}}{2}-\frac{U}{2}+\Pi+\frac{3p}{\rho^{*}}\bigg{)}d^{3}x\
$ $\displaystyle-\frac{1}{2}\int\rho^{*}W^{j}d^{3}x,$ (43)

where the post-Newtonian gravitational potential, $W^{i}$, is written as

$W^{i}=\int\frac{\rho^{*}\,[{\bf
v\cdot(x-x^{\prime})}](x-x^{\prime})^{i}}{|{\bf x}-{\bf
x^{\prime}}|^{3}}d^{3}x^{\prime}\,.$ (44)

The previous results show that 4D-EGB gravity does not violate the total
conservation of energy and momentum in the PN regime. This is an expected
outcome since the theory is Lagrangian based and, as shown in Ref. [52], they
should not violate PN conservation laws.

###  IV.1 The $\bm{\zeta}$’s and $\bm{\alpha}$’s parameters

In the context of the PPN formalism, the results obtained here directly show
that the PPN parameters $\zeta_{1},\zeta_{2},\zeta_{3},\zeta_{4}$ and
$\alpha_{3}$ do not have their physical meaning changed by the presence of the
potential $\Phi_{\mathcal{G}}$ in the PN expansion. Their values, if different
from zero, continue to indicate a violation of energy and momentum
conservation. Thus, one can conclude that, for 4D-EGB theory, one has

$\zeta_{1}=\zeta_{2}=\zeta_{3}=\zeta_{4}=\alpha_{3}=0.$ (45)

Since the expressions (36) and (43) are identical to GR counterparts, one can
also conclude that parameters $\alpha_{1}$ and $\alpha_{2}$ retain their
physical meaning unchanged since no extra-PPN potentials are present in the
expressions for the conserved energy and momentum. Moreover, a direct
comparison with the general PPN expression for $P^{j}$ assuming results in
(45), the following is obtained,

$\alpha_{1}=\alpha_{2}=0.$ (46)

The latter result shows that 4D-EGB theory does not have any preferred frame
effects (see Appendix A for details).

##  V Equation of motion of massive bodies

In this section, we split the fluid description of the source into $N$
separated bodies. We aim to obtain the PN equations of motion for the body’s
center-of-mass positions. This is a realistic way to deal with the
trajectories of massive and finite-volume bodies instead of assuming that they
can be treated as test particles. Each body indexed by $A$ has a material mass
given by

$m_{A}=\int_{A}\rho^{*}d^{3}x.$ (47)

The volume where the integration above is calculated is a time-independent
region of space that extends beyond the volume occupied by the body. It is
large enough that, in a time interval $dt$, the body does not cross its
boundary surface, but it is also small enough not to intersect with any other
body of the system. The center-of-mass of a body $A$, its velocity, and
acceleration are then defined as

$\displaystyle\bm{r}_{A}(t)\equiv$ $\displaystyle\
\dfrac{1}{m_{A}}\int\rho^{*}\bm{x}\,d^{3}x,$ (48)
$\displaystyle\bm{v}_{A}(t)\equiv$ $\displaystyle\
\dfrac{d\bm{r}_{A}}{dt}=\dfrac{1}{m_{A}}\int\rho^{*}\bm{v}\,d^{3}x,$ (49)
$\displaystyle\bm{a}_{A}(t)\equiv$ $\displaystyle\
\dfrac{d\bm{v}_{A}}{dt}=\dfrac{1}{m_{A}}\int\rho^{*}\dfrac{d\bm{v}}{dt}\,d^{3}x\,.$
(50)

The integrand of eq. (50) is found from the Euler equation with its PN
extension, which can be derived from eq. (38). Each integral is then
interpreted as a force that can be calculated according to the techniques
depicted in Ref. [46]. Ordinary GR terms give rise to 18 forces, while the
4D-EGB correction can be split into two more contributions, namely,

$\displaystyle F^{i}_{\scriptscriptstyle 19}$
$\displaystyle=\mp\,4\hat{\alpha}\int_{A}\rho^{*}\partial_{j}U\partial_{ij}U\,d^{3}x,$
(51) $\displaystyle F^{i}_{\scriptscriptstyle 20}$
$\displaystyle=\pm\,4\hat{\alpha}\int_{A}\rho^{*}\partial^{i}\Phi_{7}\,d^{3}x.$
(52)

These forces do not cancel one with each other, as in the previous section,
since the potentials in the integrand are integrals over the volume occupied
by the fluid, while the forces are integrals over the volume surrounding the
body $A$ only.

To explicitly calculate the forces above, we need to separate the potentials
into an internal part, produced by the body $A$, and an external one sourced
by the remaining bodies of the system. For instance, $U$ can be written as
follows,

$\displaystyle U$
$\displaystyle=\int_{A}\dfrac{\rho^{*}{}^{\prime}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime}+\sum_{B\neq
A}\int_{B}\dfrac{\rho^{*}{}^{\prime}}{|\bm{x}-\bm{x}^{\prime}|}\,d^{3}x^{\prime},$
$\displaystyle\equiv U_{\scriptscriptstyle A}+U_{\scriptscriptstyle
A}^{\scriptscriptstyle ext}.$ (53)

Assuming a wide separation between bodies implies that, when evaluating an
external potential within the body $A$, it can be expanded in a Taylor series.
As an example, one has

$U_{A}^{ext}(t,\bm{x})\approx
U_{A}^{ext}(t,\bm{r}_{A})+\bar{x}^{j}\partial_{j}U_{A}^{ext}(t,\bm{r}_{A})+\dots\,,$
(54)

where $\bm{\bar{x}}=\bm{x}-\bm{r}_{A}$. The series is truncated in its second
term since the next terms are at least of order ${\cal{O}}(R_{A}/r_{AB})^{2}$,
where $R_{A}$ is the typical body radius and $r_{AB}=|\bm{r}_{A}-\bm{r}_{B}|$
is the interbody distance. This expansion is used to extract the external
pieces of potentials from the integrals. Moreover, each body of the system is
assumed to be reflection-symmetric above its own center of mass, i.e.,
$\rho^{*}(t,\bm{r}_{A}-\bm{\bar{x}})=\rho^{*}(t,\bm{r}_{A}+\bm{\bar{x}})$.
This symmetry allows us to eliminate any integral having an odd number of
internal vectors, such as $\bm{\bar{x}}$.

Starting with $F^{i}_{\scriptscriptstyle 20}$, one has

$F^{i}_{\scriptscriptstyle
20}=\pm\,4\hat{\alpha}\int\rho^{*}\partial^{i}\Phi_{\scriptscriptstyle
7,A}d^{3}x\pm\,4\hat{\alpha}m_{A}\partial^{i}\Phi_{7}^{\scriptscriptstyle
ext}(t,\bm{r}_{A}).$ (55)

Using the brief notation $s=|\bm{x}-\bm{x}^{\prime}|$ one can write
$\partial_{i}\Phi_{\scriptscriptstyle 7,A}=\int_{\scriptscriptstyle
A}\rho^{*}{}^{\prime}\partial_{j}^{\prime}U^{\prime}\partial_{ij}s^{-1}\,d^{3}x^{\prime}$
and, after separating the potential inside the integrand in its internal and
external parts, one obtains,

$\displaystyle\int\rho^{*}\partial_{i}\Phi_{\scriptscriptstyle 7,A}d^{3}x=$
$\displaystyle\int_{\scriptscriptstyle
A}\rho^{*}\partial_{j}U_{\scriptscriptstyle
A}\partial_{ij}U_{\scriptscriptstyle A}d^{3}\bar{x}$ $\displaystyle\ +\
\partial_{j}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}\int_{\scriptscriptstyle A}\rho^{*}\partial_{ij}U_{\scriptscriptstyle
A}\,d^{3}\bar{x}.$ (56)

In the above expression, we have changed the variables $\bm{x}$ by
$\bm{\bar{x}}$, and used, when necessary, the reflection-symmetric properties
discussed before.

For $F^{i}_{\scriptscriptstyle 19}$, one has

$\displaystyle F^{i}_{\scriptscriptstyle 19}$
$\displaystyle=\int_{\scriptscriptstyle
A}\rho^{*}\partial_{j}U_{\scriptscriptstyle
A}\partial_{ij}U_{\scriptscriptstyle
A}d^{3}\bar{x}+\partial_{j}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}\int_{\scriptscriptstyle A}\rho^{*}\partial_{ij}U_{\scriptscriptstyle
A}\,d^{3}\bar{x}$ $\displaystyle+\partial_{ij}U^{\scriptscriptstyle
ext}_{\scriptscriptstyle A}\int_{\scriptscriptstyle
A}\rho^{*}\partial_{j}U_{\scriptscriptstyle
A}\,d^{3}\bar{x}+m_{\scriptscriptstyle A}\partial_{j}U^{\scriptscriptstyle
ext}_{\scriptscriptstyle A}\partial_{ij}U^{\scriptscriptstyle
ext}_{\scriptscriptstyle A}.$ (57)

We then note that the first line above is equal to (56) and will cancel when
both forces are added. The third integral in the r.h.s. above vanishes since
it has an odd number of internal vectors $\bar{x}$ and $\bar{x}^{\prime}$ due
to the gradient of $U$. At the end, one has,

$F^{i}_{\scriptscriptstyle 19}+F^{i}_{\scriptscriptstyle
20}=\pm\,4\hat{\alpha}m_{\scriptscriptstyle
A}\left(\partial^{i}\Phi^{\scriptscriptstyle ext}_{7,{\scriptscriptstyle
A}}-\partial_{j}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}\partial_{ij}U^{\scriptscriptstyle ext}_{\scriptscriptstyle A}\right).$ (58)

Here, the external potentials are always calculated at $\bm{x}=\bm{r}_{A}$
after the derivatives are done.

Now, we have to express the external potentials in terms of the center-of-mass
positions. We start by calculating $\partial_{j}U^{\scriptscriptstyle
ext}_{\scriptscriptstyle A}$ and $\partial_{ij}U^{\scriptscriptstyle
ext}_{\scriptscriptstyle A}$. These potentials depend on derivatives of $s$
and integrations over the volume occupied by the bodies $B\neq A$, for
instance $\partial_{j}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}=\sum_{\scriptscriptstyle B\neq
A}\int\rho^{*}{}^{\prime}\partial_{j}s^{-1}d^{3}x^{\prime}$. TThus, due to the
wide separation between bodies in the system, we can write
$s=|\bm{x}-\bm{r}_{\scriptscriptstyle B}-\bm{\bar{x}}^{\prime}|$ and expand
around $\bm{\bar{x}}=0$ in order to obtain the following derivatives,

$\displaystyle\partial_{j}s^{-1}$
$\displaystyle\approx-\frac{s^{j}_{\scriptscriptstyle
B}}{s_{\scriptscriptstyle B}^{3}}+\bar{x}^{\prime
k}\left(\frac{s^{j}_{\scriptscriptstyle B}s^{k}_{\scriptscriptstyle
B}}{s_{\scriptscriptstyle B}^{5}}-\dfrac{\delta^{jk}}{s_{\scriptscriptstyle
B}^{3}}\right),$ (59) $\displaystyle\partial_{ij}s^{-1}$
$\displaystyle\approx-\frac{\delta_{ij}}{s_{\scriptscriptstyle
B}^{3}}+3\dfrac{s^{i}_{\scriptscriptstyle B}s^{j}_{\scriptscriptstyle
B}}{s^{5}_{\scriptscriptstyle B}}$ $\displaystyle\quad+3\bar{x}^{\prime
k}\left(\dfrac{5s^{i}_{\scriptscriptstyle B}s^{j}_{\scriptscriptstyle
B}s^{k}_{\scriptscriptstyle B}}{s^{7}_{\scriptscriptstyle
B}}-\dfrac{\delta^{k(i}s^{j)}_{\scriptscriptstyle
B}+\delta^{ij}s^{k}_{\scriptscriptstyle B}}{s^{5}_{\scriptscriptstyle
B}}\right),$ (60)

where we used the short notation $\bm{s}_{\scriptscriptstyle
B}=\bm{x}-\bm{r}_{\scriptscriptstyle B}$ with $s_{\scriptscriptstyle B}$ being
its norm. Also, indices between parenthesis indicate their symmetric
counterpart. With the expansions above at hand, it is easy to verify that,

$\displaystyle\partial_{j}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}(t,\bm{r}_{\scriptscriptstyle A})$ $\displaystyle\approx-\sum_{B\neq
A}\dfrac{m_{\scriptscriptstyle B}n^{j}_{\scriptscriptstyle
AB}}{r_{\scriptscriptstyle AB}^{2}}$ (61)
$\displaystyle\partial_{ij}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}(t,\bm{r}_{\scriptscriptstyle A})$ $\displaystyle\approx-\sum_{B\neq
A}\frac{m_{\scriptscriptstyle B}}{r_{\scriptscriptstyle
AB}^{3}}\left(\delta_{ij}-3n^{i}_{\scriptscriptstyle
AB}n^{j}_{\scriptscriptstyle AB}\right),$ (62)

where $n^{i}_{\scriptscriptstyle AB}=r^{i}_{\scriptscriptstyle
AB}/r_{\scriptscriptstyle AB}$. Using the previous above results, one can
write,

$\displaystyle\partial_{j}U^{\scriptscriptstyle ext}_{\scriptscriptstyle A}$
$\displaystyle\partial_{ij}U^{\scriptscriptstyle ext}_{\scriptscriptstyle
A}=\sum_{B\neq A}\Bigg{[}-\dfrac{2m_{\scriptscriptstyle
B}^{2}}{r_{\scriptscriptstyle AB}^{5}}\,n_{\scriptscriptstyle AB}^{i}$
$\displaystyle+\\!\\!\\!\sum_{C\neq A,B}\dfrac{m_{\scriptscriptstyle
B}m_{\scriptscriptstyle C}}{r_{\scriptscriptstyle AB}^{2}r_{\scriptscriptstyle
AC}^{3}}\big{(}n_{\scriptscriptstyle AB}^{i}-3(\bm{n}_{\scriptscriptstyle
AB}\cdot\bm{n}_{\scriptscriptstyle AC})n_{\scriptscriptstyle
AC}^{i}\big{)}\Bigg{]}.$ (63)

The case of $\partial^{i}\Phi^{\scriptscriptstyle ext}_{7,{\scriptscriptstyle
A}}$ is more involved since it depends on another potential. This last must be
separated into a part generated by body $A$, another by the body $B$ (the same
as the integral limit), and a third one sourced by the remaining bodies of the
system. By doing that, one can write,

$\displaystyle\partial_{j}\Phi^{\scriptscriptstyle ext}_{7,{\scriptscriptstyle
A}}=$ $\displaystyle\sum_{B\neq A}\bigg{[}\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\partial^{\prime}_{i}U^{\prime}_{\scriptscriptstyle
A}\partial_{ij}s\,d^{3}x^{\prime}+\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\partial^{\prime}_{i}U^{\prime}_{\scriptscriptstyle
B}\partial_{ij}s\,d^{3}x^{\prime}$ $\displaystyle\quad\qquad+\sum_{C\neq
A,B}\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\partial^{\prime}_{i}U^{\prime}_{\scriptscriptstyle
C}\partial_{ij}s\,d^{3}x^{\prime}\ \bigg{]}.$ (64)

The first integral reads,

$\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\partial^{\prime}_{i}U^{\prime}_{\scriptscriptstyle
A}\partial_{ij}s\,d^{3}x=\\!\int_{\scriptscriptstyle\\!A}\\!\int_{\scriptscriptstyle\\!B}\\!\\!\rho^{*}{}^{\prime}\rho^{*}{}^{\prime\prime}\partial^{\prime}_{i}s^{\prime}\partial_{ij}s\,d^{3}x^{\prime}d^{3}x^{\prime\prime},$
(65)

with $s^{\prime}=|\bm{x}^{\prime}-\bm{x}^{\prime\prime}|$. *The expansion of
$\partial_{ij}s$ will be equivalent to (60)But the expansion of
$\partial^{\prime}_{i}s^{\prime}$ is not straightforward since
$\bm{x}^{\prime}$ is expanded around $\bm{r}_{\scriptscriptstyle B}$ whereas
$\bm{x}^{\prime\prime}$ is expanded around $\bm{r}_{\scriptscriptstyle A}$. We
only need to get the leading order of this expansion, which gives
$\partial^{\prime}_{i}s^{\prime}\approx n^{i}_{\scriptscriptstyle
AB}/r_{\scriptscriptstyle AB}^{2}$, and the first term in the r.h.s. of (64)
can be accommodated as,

$\sum_{B\neq A}\dfrac{2m_{\scriptscriptstyle A}m_{\scriptscriptstyle
B}}{r^{5}_{\scriptscriptstyle AB}}n^{j}_{\scriptscriptstyle AB}.$ (66)

The second integral in (64) is similar to (65), but now, both limits of
integration refer to the volume occupied by the body $B$. Thus, we simply
expand $\partial_{ij}s$ provided
$\partial^{\prime}_{i}s^{\prime}=(\bar{x}^{\prime}-\bar{x}^{\prime\prime})^{i}/|\bar{x}^{\prime}-\bar{x}^{\prime\prime}|^{3}$.
Using (60) once again, one arrive at the following integral,

$\displaystyle\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\partial^{\prime}_{i}U^{\prime}_{\scriptscriptstyle
B}\partial_{ij}s\,$ $\displaystyle d^{3}x^{\prime}\approx
3\left(\dfrac{5s^{i}_{\scriptscriptstyle B}s^{j}_{\scriptscriptstyle
B}s^{k}_{\scriptscriptstyle B}}{s^{7}_{\scriptscriptstyle
B}}-\dfrac{\delta^{k(i}s^{j)}_{\scriptscriptstyle
B}+\delta^{ij}s^{k}_{\scriptscriptstyle B}}{s^{5}_{\scriptscriptstyle
B}}\right)$ $\displaystyle\times\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\rho^{*}{}^{\prime\prime}\dfrac{\bar{x}^{\prime
k}(\bar{x}^{\prime}-\bar{x}^{\prime\prime})^{i}}{|\bm{\bar{x}}^{\prime}-\bm{\bar{x}}^{\prime\prime}|^{3}}d^{3}\bar{x}^{\prime}d^{3}\bar{x}^{\prime\prime}.$
(67)

The integral in the second line above can be symmetrized to acquire the
structural integral shape shown below,

$\Omega^{ik}_{\scriptscriptstyle B}=-\frac{1}{2}\int_{\scriptscriptstyle
B}\rho^{*}{}^{\prime}\rho^{*}{}^{\prime\prime}\dfrac{(\bar{x}^{\prime}-\bar{x}^{\prime\prime})^{k}(\bar{x}^{\prime}-\bar{x}^{\prime\prime})^{i}}{|\bm{\bar{x}}^{\prime}-\bm{\bar{x}}^{\prime\prime}|^{3}}d^{3}\bar{x}^{\prime}d^{3}\bar{x}^{\prime\prime}.$
(68)

It is worth noting that the trace of $\Omega^{ij}_{\scriptscriptstyle B}$
gives the gravitational energy of body $B$. At the end, the second term in
(64) will read,

$\sum_{B\neq A}\dfrac{3\Omega^{ik}_{\scriptscriptstyle
B}}{r^{4}_{\scriptscriptstyle AB}}\left(\delta_{ik}n^{j}_{\scriptscriptstyle
AB}+2\delta_{ij}n^{k}_{\scriptscriptstyle AB}-5n^{i}_{\scriptscriptstyle
AB}n^{k}_{\scriptscriptstyle AB}n^{j}_{\scriptscriptstyle AB}\right).$ (69)

The third term in (64) is calculated similar to the first one, but now
$\bm{x}^{\prime}$ is expanded around $\bm{r}_{\scriptscriptstyle B}$ while
$\bm{x}^{\prime\prime}$ is expanded around $\bm{r}_{\scriptscriptstyle C}$.
Thus, $\partial^{\prime}_{i}s^{\prime}\approx-n^{i}_{\scriptscriptstyle
BC}/r_{\scriptscriptstyle BC}^{2}$ and, using the leading part of (60) with
$\bm{x}=\bm{r}_{\scriptscriptstyle A}$, one gets,

$\sum_{B\neq A}\sum_{C\neq A,B}\dfrac{m_{\scriptscriptstyle
B}m_{\scriptscriptstyle C}}{r_{\scriptscriptstyle AB}^{3}r_{\scriptscriptstyle
BC}^{2}}\left[n_{\scriptscriptstyle BC}^{j}-3(\bm{n}_{\scriptscriptstyle
AB}\cdot\bm{n}_{\scriptscriptstyle BC})n_{\scriptscriptstyle AB}^{j}\right].$
(70)

By collecting all the previous results, the external contribution to the
gradient of $\Phi^{\scriptscriptstyle ext}_{\scriptscriptstyle 7,A}$ can be
recast as

$\partial_{j}\Phi^{\scriptscriptstyle ext}_{\scriptscriptstyle
7,A}=\sum_{B\neq A}\Bigg{[}\dfrac{2m_{\scriptscriptstyle
A}m_{\scriptscriptstyle B}}{r^{5}_{\scriptscriptstyle
AB}}n^{j}_{\scriptscriptstyle AB}+\dfrac{3\Omega^{ik}_{\scriptscriptstyle
B}}{r^{4}_{\scriptscriptstyle AB}}\left(\delta_{ik}n^{j}_{\scriptscriptstyle
AB}+2\delta_{ij}n^{k}_{\scriptscriptstyle AB}-5n^{i}_{\scriptscriptstyle
AB}n^{k}_{\scriptscriptstyle AB}n^{j}_{\scriptscriptstyle
AB}\right)+\\!\\!\\!\sum_{C\neq A,B}\dfrac{m_{\scriptscriptstyle
B}m_{\scriptscriptstyle C}}{r_{\scriptscriptstyle AB}^{3}r_{\scriptscriptstyle
BC}^{2}}\big{(}n_{\scriptscriptstyle BC}^{j}-3(\bm{n}_{\scriptscriptstyle
AB}\cdot\bm{n}_{\scriptscriptstyle BC})n_{\scriptscriptstyle
AB}^{j}\big{)}\Bigg{]}.\\!$ (71)

The 4D-EGB correction to the PN-equation of motion of well-separated massive
bodies is given by the force expressions (58). Substituting in that the
results given in (V) and (71), one can explicitly write the acceleration of
body $A$ in a system of $N$ bodies as follows,

$\displaystyle a_{A}^{j}=$ $\displaystyle\ a_{\scriptscriptstyle
A}^{j}{\rm[GR]}\mp 8\hat{\alpha}\sum_{\scriptscriptstyle B\neq
A}\frac{M_{\scriptscriptstyle B}(M_{\scriptscriptstyle
A}+M_{\scriptscriptstyle A})}{r_{\scriptscriptstyle
AB}^{5}}\,n_{\scriptscriptstyle AB}^{j}\pm
6\hat{\alpha}\sum_{\scriptscriptstyle B\neq A}\frac{\Omega_{\scriptscriptstyle
B}^{ik}}{r_{\scriptscriptstyle
AB}^{4}}\left(\delta^{ik}n^{j}_{\scriptscriptstyle
AB}+2\delta^{ij}n^{k}_{\scriptscriptstyle AB}-5n^{i}_{\scriptscriptstyle
AB}n^{k}_{\scriptscriptstyle AB}n^{j}_{\scriptscriptstyle AB}\right)$
$\displaystyle\ \pm 4\hat{\alpha}\sum_{\scriptscriptstyle B\neq
A}\sum_{\scriptscriptstyle C\neq A,B}\frac{M_{\scriptscriptstyle
B}M_{\scriptscriptstyle C}}{r_{\scriptscriptstyle
AB}^{2}}\Bigg{[}\frac{n_{\scriptscriptstyle AB}^{j}}{r_{\scriptscriptstyle
AC}^{3}}-3\frac{(\bm{n}_{\scriptscriptstyle AB}\cdot\bm{n}_{\scriptscriptstyle
AC})n_{\scriptscriptstyle AC}^{j}}{r_{\scriptscriptstyle
AC}^{3}}-\frac{n_{\scriptscriptstyle BC}^{j}}{r_{\scriptscriptstyle
AB}r_{\scriptscriptstyle BC}^{2}}+3\frac{(\bm{n}_{\scriptscriptstyle
AB}\cdot\bm{n}_{\scriptscriptstyle BC})n_{\scriptscriptstyle
AB}^{j}}{r_{\scriptscriptstyle AB}r_{\scriptscriptstyle BC}^{2}}\Bigg{]}.$
(72)

In Eq. (72), $a_{\scriptscriptstyle A}^{j}{\rm[GR]}$ represents the standard
GR term to the equation of motion of each body, say body $A$, where each mass
term $m$ was replaced by $M$ provided to mass-energy is conserved [cf. (37)],
and the corrections introduced by these exchanges are beyond the PN order. In
addition to that, the terms proportional to the re-scaled coupling
$\hat{\alpha}$ accommodate the new contributions to the acceleration of the
$A$ body. In the next section, we will employ (72) to perform some
astrophysical constraints on the 4D-EGB theory, for instance, by obtaining the
orbital periastron shift in a binary system.

##  VI Astrophysical constraints

###  VI.1 The periastron advance rate

To analyze the 4D-EGB contributions to the orbital periastron advance, we
consider the two-body problem. We work in terms of relative acceleration,
$\vec{a}=\vec{a}_{1}-\vec{a}_{2}$, and use the notation
$\vec{r}=\vec{r}_{1}-\vec{r}_{2}$, $\hat{n}=\vec{r}/r$, and $m=M_{1}+M_{2}$.
We focus on the case of nearly spherical bodies; the approximation
$\Omega_{\scriptscriptstyle
N}^{ik}\approx\delta^{ik}\Omega_{\scriptscriptstyle N}/3$ remains valid.
Replacing the latter approximation in (72), one then obtains a simplified
version for the total acceleration,

$\vec{a}=\vec{a}{\rm[GR]}\pm\frac{8\hat{\alpha}m^{2}}{r^{5}}\,\hat{n}.$ (73)

The 4D-EGB contribution to the two-body system relative acceleration is of PN
order; thus, it can be considered a perturbing acceleration that gives origin
to perturbed Keplerian orbits.

Exploiting the method of osculating orbital elements (see Ref. [46]), the
perturbed motion is interpreted as an instantaneous Keplerian orbit, and the
system is always describing an ellipse, but with the orbital elements varying
with time. The perturbing acceleration is decomposed in a radial component
${\cal R}$, a component ${\cal W}$ normal to the orbital plane, and a third
component ${\cal S}$ transversal to the previous two directions. Equation (73)
thus shows,

${\cal R}=\frac{8\hat{\alpha}m^{2}}{r^{5}},\quad{\cal W}={\cal S}=0.$ (74)

The variations of the orbital elements can then be directly related to these
components.

###  VI.2 Mercury’s periastron advance rate

Once observations are made in geocentric coordinates, the periastron argument
measured is relative to the equinox. Consequently, it can be recast as
$\tilde{\omega}=\omega+\Omega\cos i$, where $\Omega$ is the angle from the
ascending node to the Earth-Sun direction, $\omega$ is the angle between the
periastron and the ascending node. At the same time, $i$ indicates the orbital
plane inclination angle relative to the ecliptic [43]. Using the fact that,
for all planets in the solar system, $i$ is too small, the variation of
$\tilde{\omega}$ then is determined by the following derivative,

$\frac{d\tilde{\omega}}{df}=\frac{p^{2}\sin f}{em(1+e\cos f)^{2}}\bigg{[}{\cal
S}\,\frac{2+e\cos f}{1+e\cos f}-{\cal R}\cot f\bigg{]},$ (75)

where $f$ is the true anomaly (the angular position of the planet with respect
to its periastron direction), $p=a(1-e^{2})$ is the semi-latus rectum, whereas
$a$ denotes the semi-major axis, and $e$ stands for the eccentricity of the
orbit. Within this Keplerian description, the Sun-planet relative position is
given by the usual conic equation,

$r=\frac{p}{1+e\cos f}.$ (76)

Integrating over a complete orbital period gives us the secular change in the
periastron position. Using (73), one can calculate the 4D-EGB contribution to
this effect, and, after including the GR contribution, one obtains

$\Delta\tilde{\omega}=\frac{6\pi
m}{p}\left[1\pm\frac{\hat{\alpha}}{p^{2}}\,(4+e^{2})\right].$ (77)

The result above is precisely found in Ref. [38], where the authors considered
the planets as described by test particles. Although the general equation of
motion (73) shows terms depending explicitly on the body’s internal structure,
through their gravitational energy tensor $\Omega^{ik}$, all those terms
cancel each other when assuming nearly spherical bodies. Therefore, it remains
to be analyzed in situations where the internal structure of the bodies does
contribute to the orbital motion.

The advance per orbit $\Delta\tilde{\omega}$ can be converted to a rate by
dividing it by the orbital period, $P$. The result is given by

$\dot{\tilde{\omega}}=3\left(\frac{2\pi}{P}\right)^{5/3}\frac{m^{2/3}}{1-e^{2}}\
\pm\
3\hat{\alpha}\left(\frac{2\pi}{P}\right)^{3}\frac{4+e^{2}}{(1-e^{2})^{2}},$
(78)

after using Kepler’s third law, $P^{2}=4\pi^{2}p^{3}/m(1-e^{2})^{3}$.

Eq. (78) can be used to put some constraints over the parameter
$\hat{\alpha}$. Data analysis from the MESSENGER mission estimates secular
Mercury’s periastron precession as $\dot{\tilde{\omega}}=(42.9799\pm 0.0009)$
arcsec per century [53]. With the help of the orbital parameters $e=0.2056$,
$P=87.97$ days, $M_{\odot}=1.9891\times 10^{30}$ kg and
$M_{\text{\mercury}}=3.3011\times 10^{23}$ kg [54], one can find an upper
bound limit for the GB coupling, which is given by

$|\hat{\alpha}|\lesssim 1.67\times 10^{16}\,\text{m}^{2}.$ (79)

This constraint is two orders of magnitude stronger than the one established
in Ref. [38] due to the use of less precise data on Mercury’s orbit.

###  VI.3 The Double Pulsar periastron shift

The binary’s orbital motion affects the pulses from a pulsar in a binary
system, altering their arrival times at radio telescopes [55]. By measuring
these arrival times accurately over long periods of time, we can detect slight
variations in the orbital motion. In 2003, the double pulsar PSR J0737 –
3039A/B was first detected, representing probably one of the best binary
systems for studying gravitational effects. Such a claim is based on the
following facts: $i)$ it is the only binary where we can see both pulsars,
$ii)$ it is not so far from Earth — nearly $0.61\text{kpc}$, and $iii)$ its
orbit is almost “edge-on” to us. All these listed properties indicate why the
double pulsar offers a unique opportunity to explore the strong gravity regime
and, consequently, any deviation from GR.

More recently Kramer  _et al._ [56] presented precise determination of
periastron precession using 16-year data span of the PSR J0737 – 3039A/B. To
obtain stringent bounds using this higher accurate system, it is necessary to
derive a new expression for $\dot{\tilde{\omega}}$, since the inclination $i$
is now close to $90^{\circ}$. Being $\cos i\approx 0$, the variations of
$\Omega$ will not affect the periastron advance rate, which, in turn, will
lead to a modification in equation (75) proportional to ${\cal W}$. However,
due to the radial nature of the 4D-EGB correction to the relative
acceleration, the expression (78) for $\dot{\tilde{\omega}}$ remains valid for
the double pulsar system. The total mass of the two stars is $2.587M_{\odot}$,
their orbit is slightly elliptical with $e=0.0878$, and they complete one
revolution in $P=0.10225$ day. The measured precession rate reads
$\dot{\tilde{\omega}}=(16.899323\pm 0.000013)\,$ degrees per year. Using these
data, it is possible to obtain a stronger constraint on the re-scaled GB
coupling, namely,

$|\hat{\alpha}|\lesssim 1.47\times 10^{11}\,\text{m}^{2}.$ (80)

By comparing the bounds found in Sections VIB and VIC, it is fair to state
that the double pulsar system imposes a much stronger constraint by reducing
the bounds on $\hat{\alpha}$ by five orders of magnitude. Nevertheless, the
above estimation on GB coupling represents an interesting improvements in
regard to the previous best estimates, $|\hat{\alpha}|<10^{15}\,\text{m}^{2}$
[38].

###  VI.4 Other observational bounds

Let’s discuss the improvement introduced on the $\hat{\alpha}$ coupling’s
bound concerning other observational tests reported in the literature in
different regimes and sectors of the 4D-EGB theory. Even though the
constraints coming from the double pulsar periastron shift enhanced the bounds
on $\hat{\alpha}$ by several orders of magnitude, it fell short regarding the
bounds emerging in the strong-field regime provided by the binary black hole
system [38]. In other words,
$\hat{\alpha}_{BBH}<\hat{\alpha}_{DPP}<10^{10}\,{\rm{m}}^{2}$. The latter
situation is considerably improved once is assumed that lightest component of
the GW190814 event could potentially correspond to stellar black hole with
mass $M=2.59^{+0.08}_{-0.09}M_{\odot}$, yielding
$\hat{\alpha}_{GW190814}<59\times 10^{6}\,{\rm{m}}^{2}$ [40]. When we account
for an extra-dark radiation term at both the background and perturbative
levels in cosmology, we can constrain it using the ACTPol data alone. This
gives us $\hat{\alpha}C^{2}=(-9\pm 6)\times 10^{-6}H^{2}_{0}$, where
$\hat{\alpha}C^{2}=\Omega^{\rm{extra}}_{\rm{rad}}H^{2}_{0}$ [42]. If we ignore
the geometrical dark radiation term, however, we get a less robust estimate on
$\hat{\alpha}<10^{36}\,{\rm{m}}^{2}$ [37]. Given that nucleosynthesis started
at an energy scale of $E\sim\text{Mev}$ when the universe is dominated by
radiation, we can estimate that $\hat{\alpha}<10^{18}\,{\rm{m}}^{2}$ [37].
However, this bound is much weaker than the double pulsar’s best estimates.

##  VII The $\bm{\beta}$ parameter and the new EGB parameter: an extended PPN
version

The previous sections showed that 4D-EGB corrections to GR post-Newtonian
metric, does not change how light travels or how energy and momentum are
conserved in the weak gravity limit up to the first PN order. Within the PPN
formalism, these results demonstrate that the physical meaning associated with
parameters $\gamma$, $\zeta$’s, and $\alpha$’s are not influenced by the new
potential $\Phi_{\mathcal{G}}$. Consequently, once it was shown that photon
geodesics and total energy and momentum expressions are the same as in GR, it
was possible to conclude that $\gamma=1$ and all $\zeta$’s and $\alpha$’s are
null. The remaining PPN parameters to be determined are $\beta$ and $\xi$. To
move forward with the determination of the PPN parameters, one assumes the
knowledge of the already determined parameters in juxtaposition with the
general PPN expression for $\Delta\omega$, namely,

$\displaystyle\Delta\tilde{\omega}_{\rm PPN}=\frac{2\pi m}{p}(4-\beta),$ (81)

one might be tempting, by comparison with (77), to define an effective $\beta$
parameter that would encompass the 4D-EGB corrections to GR. However, we
understand that such treatment is not appropriate once it would lead to a
system-dependent parameter, a rather strange feature for the PPN formalism.
More drastically, if an effective parameter is defined as a system-dependent
constant, one could not simply substitute the original parameter with its
effective version in each observable physical phenomenon. As an example, we
emphasize that if an effective $\beta$ were defined through (81) and (77),
this new parameter would not have any influence on the so-called Nordtvedt
effect: the violation of the weak equivalence principle due to explicit
contributions of self-gravitational energy to a body’s inertial and
gravitational mass.

The Nordtvedt effect can be tested for the Earth-Moon system by studying its
motion in the Sun’s gravitational field. If there is any difference between
the Earth’s free-fall acceleration and the Moon’s towards the Sun, this effect
will be parametrized by the Nordtvedt parameter,

$\eta=4\beta-4-\frac{10}{3}\,\xi,$ (82)

where $\xi$ is the PPN parameter related to the existence of preferred-
locations effects, and, once again, we use $\gamma=1$ and
$\zeta_{i}=\alpha_{i}=0$. It has been shown that $\Phi_{\mathcal{G}}$ does not
bring any dependence in the acceleration of a body $A$ with its internal
structure up to the Newtonian order [cf. Eq. (50)]. Consequently,$\eta=0$ and
no Nordtvedt effect is present, just like in GR. Therefore, a simple
substitution of $\beta$, in (82), by an effective parameter obtained from the
periastron advance would lead to erroneous conclusions for any theory.

The correct approach to tackle the issue above would be to propose a new
parameter, say $\tilde{\beta}$, to quantify the dynamical effects associated
with $\Phi_{\mathcal{G}}$. In this extended PPN version (EPPN), the Nordtvedt
parameter will remain the same, while the periastron advance per orbit will be
given by

$\displaystyle\Delta\tilde{\omega}_{\rm EPPN}=\Delta\tilde{\omega}_{\rm
PPN}+\frac{\tilde{\beta}}{p^{2}}\,(4+e^{2}).$ (83)

Within general relativity, for instance, $\beta=1$ and $\tilde{\beta}=0$,
while in the case of 4D-EGB theory one has $\beta=1$ and
$\tilde{\beta}=\pm\,\hat{\alpha}$. Hence, with $\eta=0$ one would also obtain
$\xi=0$. However, it is important to emphasize that, within this extended PPN
approach, $\beta$ alone does not determine the post-Newtonian periastron
advance effect. For instance, the $\beta=1$ result is insufficient to ensure
agreement with Mercury’s periastron advance test. Still, it is essential to
guarantee no violations of the weak equivalence principle in the case of fully
conservative theories.

##  VIII Summary

We presented a comprehensive analysis of the post-Newtonian version of the
regularized four-dimensional Einstein-Gauss-Bonnet gravitational theory. Our
investigation involved expanding the metric and the extra scalar field around
the Minkowski background and the constant scalar field value, respectively. To
model the matter content, we adopted a perfect fluid approach. The post-
Newtonian metric components are detailed in equations (14) through (16).
Notably, the primary correction emerges at fourth order in the $g_{00}$
component. Unlike previous works within the post-Newtonian formalism, we
utilized the harmonic gauge and expressed the gravitational potential in terms
of the conserved density.

Starting from conservation laws, it was shown that there is not violation of
energy and/or momentum conservation within the PN limit, confirming that the
four of $\zeta_{i}$ PPN parameters and $\alpha_{3}$ vanish. Moreover, the
expressions for energy and momentum showed that parameters $\alpha_{1}$ and
$\alpha_{2}$ are also zero; this highlights that the regularized four-
dimensional Einstein-Gauss-Bonnet gravitational theory lacks any preferred
frame effects.

We investigated the equation of motion for massive bodies and derived the
general expression for a body within a system of N bodies. Furthermore, we
recovered the standard GR term and identified several contributions arising
from the Gauss-Bonnet (GB) term. The latter discovery involves lengthy
computations, but it was essential for imposing physical bounds on the re-
scaled GB coupling.

From an observational perspective, we examined the two-body problem for nearly
spherical objects and utilized the method of osculating orbital elements to
derive the PN corrections in the periastron advance rate due to the Gauss-
Bonnet term. By analyzing data from the MESSENGER mission applied to secular
Mercury’s periastron precession, we arrived at an updated estimate for the GB
coupling: $|\hat{\alpha}|\lesssim 1.67\times 10^{16}\,\text{m}^{2}$;
consistent with previous results [38]. Further, we focused on the double
pulsar J0737 – 3039A/B, comprising two active radio pulsars in a binary system
with a period of 2.45 hours and a mild eccentricity of $e=0.088$. Taking into
account the data collected on extensive pulse timing experiments over 16 years
on six different radio telescopes around the world, this experiment reported
the precession rate for the double pulsar system,
$\dot{\tilde{\omega}}=(16.899323\pm 0.000013)\,$ degrees per year (see Ref.
[56] for further details). We proceeded by using this precession rate of the
double pulsar system to put a tighter constraint on the GB coupling, yielding
$|\hat{\alpha}|\lesssim 1.47\times 10^{11}\,\text{m}^{2}$; reducing the former
bound in five orders of magnitude. We also compared other astrophysical and
cosmological bounds on $\hat{\alpha}$ to contextualize the newly updated
results reported here.

Finally, we discussed, to a certain degree, the $\beta$ PPN parameter
identification and the introduction of a new EGB parameter. This analysis
focused on the different consequences of misreading the $\beta$ PPN parameter
within a modified gravity scenario such as the 4D-EGB model.

##  IX Acknowledgments

M.G.R. was partially supported by FAPES (Brazil), CNPq (Brazil), and by Wölfe
Foundation.

##  Appendix A The PPN formalism – a synthesis

In order to distinguish between alternative theories of gravity within the
post-Newtonian regime, the PPN formalism assumes a generic metric expansion
under reasonable (but restrictive) assumptions over the possible gravitational
potentials to be present.

The original PPN metric, as presented in Ref. [43], is written in the so-
called standard gauge. Nowadays, it is more usual for PN analysis to work with
the harmonic gauge. Although the gauge choice is a convenience matter, the
harmonic gauge usage allows to consistently connect the post-Newtonian
approximation to eh wider scopus of a post-Minkowskian theory (see Ref. [46]
for details). In the harmonic gauge, the PPN metric reads,

$\displaystyle g_{00}=-$ $\displaystyle 1+2U-2\beta
U^{2}+(2\gamma+2+\alpha_{3}+\zeta_{1}-2\xi)\Phi_{1}\,+$ $\displaystyle+$
$\displaystyle\,2(3\gamma-2\beta+1+\zeta_{2}+\xi)\Phi_{2}\
+2(1+\zeta_{3})\Phi_{3}\,+$ $\displaystyle+$
$\displaystyle\,2(3\gamma+3\zeta_{4}-2\xi)\Phi_{4}-(\zeta_{1}-2\xi)\Phi_{6}\,-$
$\displaystyle-$
$\displaystyle\,2\xi\Phi_{W}+(1+\alpha_{2}-\zeta_{1}+2\xi)\partial_{tt}\chi$
$\displaystyle g_{0i}=-$ $\displaystyle(2\gamma+2+\alpha_{1}/2)V_{i}$ (84)
$\displaystyle g_{ij}=\ ($ $\displaystyle 1+2\gamma\,U)\,\delta_{ij}\,.$

The organization of the PPN parameters in the metric above is done in such a
way that is possible to give to each one of them a specific physical meaning
which associates them to a property or a measurable effect. Table 1 the
significance of each PPN parameter and the current bounds on PPN parameters
are shown in Table 2.

Table 1: The PPN Parameters physical meaning. Parameter |  Physical meaning associated | Value in GR
---|---|---
$\gamma$ |  light motion tests | 1
$\beta$ |  Periastron shift in a binary system | 1
$\xi$ |  Preferred-location effects | 0
$\alpha_{1}$, $\alpha_{2}$ |  Preferred-frame effects | all null
and $\alpha_{3}$ |  |
$\zeta_{1}$, $\zeta_{2}$, $\zeta_{3}$, |  Violation of conservation | all null
$\zeta_{4}$ and $\alpha_{3}$ |  of total momentum |
Table 2: Current limits on PPN parameters (from Ref. [1]).  Parameter | LimitAA | Parameter | Limit
---|---|---|---
$\gamma-1$ |  ${\scriptstyle\lesssim\,}2.3\times 10^{-5\phantom{0}}$ | $\alpha_{3}$ |  ${\scriptstyle\lesssim\,}4\times 10^{-20}$
$\beta-1$ |  ${\scriptstyle\lesssim\,}8\times 10^{-5\phantom{0}}$ | $\zeta_{1}$ |  ${\scriptstyle\lesssim\,}2\times 10^{-2\phantom{0}}$
$\xi$ |  ${\scriptstyle\lesssim\,}4\times 10^{-9\phantom{0}}$ | $\zeta_{2}$ |  ${\scriptstyle\lesssim\,}4\times 10^{-5\phantom{0}}$
$\alpha_{1}$ |  ${\scriptstyle\lesssim\,}7\times 10^{-5\phantom{0}}$ | $\zeta_{3}$ |  ${\scriptstyle\lesssim\,}10^{-8\phantom{0}}$
$\alpha_{2}$ |  ${\scriptstyle\lesssim\,}2\times 10^{-9\phantom{0}}$ | $\zeta_{4}$ | —

When a theory does not present violations of energy and momentum at the PN
regime, all $\zeta$’s and $\alpha_{3}$ vanish. The expressions for the
conserved energy and momentum are then reduced to

$E\equiv\int\left(\frac{1}{2}\rho^{*}v^{2}+\rho^{*}\Pi-\frac{1}{2}\rho^{*}U\right)d^{3}x.$
(85) $\displaystyle P_{j}=$
$\displaystyle\int\rho^{\ast}v_{j}\left(1+\frac{1}{2}v^{2}-\frac{1}{2}U+\Pi+\frac{p}{\rho^{\ast}}\right)d^{3}x$
(86)
$\displaystyle-\frac{1}{2}\int\rho^{\ast}\left[(1+\alpha_{2})\Phi_{j}+\frac{1}{2}(\alpha_{1}-\alpha_{2})V_{j}\right]d^{3}x,$

## References

  * Will [2014] C. M. Will, The Confrontation between General Relativity and Experiment, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Ishak [2019] M. Ishak, Testing General Relativity in Cosmology, Living Reviews in Relativity 22, 1 (2019), arXiv:1806.10122 [astro-ph.CO] .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Riess  _et al._ [1998] A. G. Riess _et al._ (Supernova Search Team), Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ [1999] S. Perlmutter _et al._ (Supernova Cosmology Project), Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae, Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Weinberg [1989] S. Weinberg, The Cosmological Constant Problem, Rev. Mod. Phys. 61, 1 (1989).
  * Carroll [2001] S. M. Carroll, The Cosmological constant, Living Rev. Rel. 4, 1 (2001), arXiv:astro-ph/0004075 .
  * Peebles and Ratra [2003] P. J. E. Peebles and B. Ratra, The Cosmological Constant and Dark Energy, Rev. Mod. Phys. 75, 559 (2003), arXiv:astro-ph/0207347 .
  * Frieman  _et al._ [2008] J. Frieman, M. Turner, and D. Huterer, Dark Energy and the Accelerating Universe, Ann. Rev. Astron. Astrophys. 46, 385 (2008), arXiv:0803.0982 [astro-ph] .
  * Martin [2008] J. Martin, Quintessence: a mini-review, Mod. Phys. Lett. A 23, 1252 (2008), arXiv:0803.4076 [astro-ph] .
  * Tsujikawa [2013] S. Tsujikawa, Quintessence: A Review, Class. Quant. Grav. 30, 214003 (2013), arXiv:1304.1961 [gr-qc] .
  * Armendariz-Picon  _et al._ [2000] C. Armendariz-Picon, V. F. Mukhanov, and P. J. Steinhardt, A Dynamical solution to the problem of a small cosmological constant and late time cosmic acceleration, Phys. Rev. Lett. 85, 4438 (2000), arXiv:astro-ph/0004134 .
  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified Gravity and Cosmology, Phys.Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Joyce  _et al._ [2016] A. Joyce, L. Lombriser, and F. Schmidt, Dark Energy Versus Modified Gravity, Ann. Rev. Nucl. Part. Sci. 66, 95 (2016), arXiv:1601.06133 [astro-ph.CO] .
  * Nojiri  _et al._ [2017] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Heisenberg [2019] L. Heisenberg, A systematic approach to generalisations of General Relativity and their cosmological implications, Phys. Rept. 796, 1 (2019), arXiv:1807.01725 [gr-qc] .
  * Lovelock [1972] D. Lovelock, The four-dimensionality of space and the einstein tensor, J. Math. Phys. 13, 874 (1972).
  * Langlois and Noui [2016] D. Langlois and K. Noui, Degenerate higher derivative theories beyond Horndeski: evading the Ostrogradski instability, JCAP 02, 034, arXiv:1510.06930 [gr-qc] .
  * de Rham and Matas [2016] C. de Rham and A. Matas, Ostrogradsky in Theories with Multiple Fields, JCAP 06, 041, arXiv:1604.08638 [hep-th] .
  * Horndeski [1974] G. W. Horndeski, Second-order scalar-tensor field equations in a four-dimensional space, Int.J.Theor.Phys. 10, 363 (1974).
  * Gleyzes  _et al._ [2015] J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, Healthy theories beyond Horndeski, Phys. Rev. Lett. 114, 211101 (2015), arXiv:1404.6495 [hep-th] .
  * Barausse and Yagi [2015] E. Barausse and K. Yagi, Gravitation-Wave Emission in Shift-Symmetric Horndeski Theories, Phys. Rev. Lett. 115, 211105 (2015), arXiv:1509.04539 [gr-qc] .
  * Kase and Tsujikawa [2019] R. Kase and S. Tsujikawa, Dark energy in Horndeski theories after GW170817: A review, Int. J. Mod. Phys. D 28, 1942005 (2019), arXiv:1809.08735 [gr-qc] .
  * Kobayashi [2019] T. Kobayashi, Horndeski theory and beyond: a review, Rept. Prog. Phys. 82, 086901 (2019), arXiv:1901.07183 [gr-qc] .
  * Ezquiaga and Zumalacárregui [2017] J. M. Ezquiaga and M. Zumalacárregui, Dark Energy After GW170817: Dead Ends and the Road Ahead, Phys. Rev. Lett. 119, 251304 (2017), arXiv:1710.05901 [astro-ph.CO] .
  * Creminelli and Vernizzi [2017] P. Creminelli and F. Vernizzi, Dark Energy after GW170817 and GRB170817A, Phys. Rev. Lett. 119, 251302 (2017), arXiv:1710.05877 [astro-ph.CO] .
  * Fernandes  _et al._ [2020] P. G. Fernandes, P. Carrilho, T. Clifton, and D. J. Mulryne, Derivation of Regularized Field Equations for the Einstein-Gauss-Bonnet Theory in Four Dimensions, Phys. Rev. D 102, 024025 (2020), arXiv:2004.08362 [gr-qc] .
  * Hennigar  _et al._ [2020] R. A. Hennigar, D. Kubizňák, R. B. Mann, and C. Pollack, On taking the D → 4 limit of Gauss-Bonnet gravity: theory and solutions, JHEP 07, 027, arXiv:2004.09472 [gr-qc] .
  * Lu and Pang [2020] H. Lu and Y. Pang, Horndeski gravity as $D\rightarrow 4$ limit of Gauss-Bonnet, Phys. Lett. B 809, 135717 (2020), arXiv:2003.11552 [gr-qc] .
  * Kobayashi [2020] T. Kobayashi, Effective scalar-tensor description of regularized Lovelock gravity in four dimensions, JCAP 07, 013, arXiv:2003.12771 [gr-qc] .
  * Glavan and Lin [2020] D. Glavan and C. Lin, Einstein-Gauss-Bonnet Gravity in Four-Dimensional Spacetime, Phys. Rev. Lett. 124, 081301 (2020), arXiv:1905.03601 [gr-qc] .
  * Gurses  _et al._ [2020] M. Gurses, T. c. Şişman, and B. Tekin, Comment on ”Einstein-Gauss-Bonnet Gravity in 4-Dimensional Space-Time”, Phys. Rev. Lett. 125, 149001 (2020), arXiv:2009.13508 [gr-qc] .
  * Arrechea  _et al._ [2020] J. Arrechea, A. Delhom, and A. Jiménez-Cano, Comment on “Einstein-Gauss-Bonnet Gravity in Four-Dimensional Spacetime”, Phys. Rev. Lett. 125, 149002 (2020), arXiv:2009.10715 [gr-qc] .
  * Ai [2020] W.-Y. Ai, A note on the novel 4D Einstein–Gauss–Bonnet gravity, Commun. Theor. Phys. 72, 095402 (2020), arXiv:2004.02858 [gr-qc] .
  * Fernandes  _et al._ [2022] P. G. S. Fernandes, P. Carrilho, T. Clifton, and D. J. Mulryne, The 4D Einstein–Gauss–Bonnet theory of gravity: a review, Class. Quant. Grav. 39, 063001 (2022), arXiv:2202.13908 [gr-qc] .
  * Clifton  _et al._ [2020] T. Clifton, P. Carrilho, P. G. S. Fernandes, and D. J. Mulryne, Observational constraints on the regularized 4D Einstein-Gauss-Bonnet theory of gravity, Phys. Rev. D 102, 084005 (2020).
  * Feng  _et al._ [2021] J.-X. Feng, B.-M. Gu, and F.-W. Shu, Theoretical and observational constraints on regularized 4$D$ Einstein-Gauss-Bonnet gravity, Phys. Rev. D 103, 064002 (2021), arXiv:2006.16751 [gr-qc] .
  * Charmousis  _et al._ [2022] C. Charmousis, A. Lehébel, E. Smyrniotis, and N. Stergioulas, Astrophysical constraints on compact objects in 4D Einstein-Gauss-Bonnet gravity, JCAP 02 (02), 033, arXiv:2109.01149 [gr-qc] .
  * Aoki  _et al._ [2020] K. Aoki, M. A. Gorji, and S. Mukohyama, Cosmology and gravitational waves in consistent $D\to 4$ Einstein-Gauss-Bonnet gravity, JCAP 09, 014, [Erratum: JCAP 05, E01 (2021)], arXiv:2005.08428 [gr-qc] .
  * Zanoletti  _et al._ [2024] C. M. A. Zanoletti, B. R. Hull, C. D. Leonard, and R. B. Mann, Cosmological constraints on 4-dimensional Einstein-Gauss-Bonnet gravity, JCAP 01, 043, arXiv:2310.19871 [astro-ph.CO] .
  * Will [1993] C. M. Will,  _Theory and experiment in gravitational physics_ , 2nd ed. (Cambridge University Press, Cambridge, 1993).
  * Easson  _et al._ [2020] D. A. Easson, T. Manton, and A. Svesko, $D\to 4$ Einstein-Gauss-Bonnet gravity and beyond, JCAP 10, 026, arXiv:2005.12292 [hep-th] .
  * Hohmann [2015] M. Hohmann, Parametrized post-Newtonian limit of Horndeski’s gravity theory, Phys. Rev. D92, 064019 (2015), arXiv:1506.04253 [gr-qc] .
  * Poisson and Will [2014] E. Poisson and C. M. Will,  _Gravity: Newtonian, Post-Newtonian, Relativistic_ (Cambridge University Press, 2014).
  * Bittencourt  _et al._ [2014] E. Bittencourt, U. Moschella, M. Novello, and J. D. Toniato, Cosmology in geometric scalar gravity, Phys. Rev. D90, 123540 (2014), arXiv:1412.4227 [gr-qc] .
  * Toniato  _et al._ [2017] J. D. Toniato, D. C. Rodrigues, A. O. F. de Almeida, and N. Bertini, Will-Nordtvedt PPN formalism applied to renormalization group extensions of general relativity, Phys. Rev. D96, 064034 (2017), arXiv:1706.09032 [gr-qc] .
  * Toniato  _et al._ [2020] J. D. Toniato, D. C. Rodrigues, and A. Wojnar, Palatini $f(R)$ gravity in the solar system: post-Newtonian equations of motion and complete PPN parameters, Phys. Rev. D 101, 064050 (2020), arXiv:1912.12234 [gr-qc] .
  * Toniato and Rodrigues [2021] J. D. Toniato and D. C. Rodrigues, Post-Newtonian $\gamma$-like parameters and the gravitational slip in scalar-tensor and f(R) theories, Phys. Rev. D 104, 044020 (2021), arXiv:2106.12542 [gr-qc] .
  * Alves  _et al._ [2024] M. F. S. Alves, J. D. Toniato, and D. C. Rodrigues, A detailed first-order post-newtonian analysis of massive brans-dicke theories: numerical constraints and the $\beta$ parameter meaning (2024), accepted for publication in Physical Review D, arXiv:2307.11883 [gr-qc] .
  * Lee  _et al._ [1974] D. L. Lee, A. P. Lightman, and W. T. Ni, Conservation laws and variational principles in metric theories of gravity, Phys. Rev. D10, 1685 (1974).
  * Park  _et al._ [2017] R. S. Park, W. M. Folkner, A. S. Konopliv, J. G. Williams, D. E. Smith, and M. T. Zuber, Precession of mercury’s perihelion from ranging to the messenger spacecraft, The Astronomical Journal 153, 121 (2017).
  * Stark  _et al._ [2015] A. Stark, J. Oberst, and H. Hussmann, Mercury’s resonant rotation from secular orbital elements, Celestial Mechanics and Dynamical Astronomy 123, 263 (2015).
  * Taylor [1994] J. H. Taylor, Binary pulsars and relativistic gravity, Rev. Mod. Phys. 66, 711 (1994).
  * Kramer  _et al._ [2021] 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).