# The Equivalence and/or the Effacing principle in $f\left(R\right)$ theories
of gravity

Soham Bhattacharyya  Max Planck Institute for Gravitational Physics (Albert
Einstein Institute)
Leibniz Universität Hannover<EMAIL_ADDRESS>
###### Abstract

The Einstein-Hilbert action of general theory of relativity (GR) is the
integral of the scalar curvature $R$. It is a theory that is drawn from the
Equivalence principle, and has predictions that come out as a consequence of
the principle, in observables. Testing such observables to find
confirmation/infirmation of the principle have formed a significant chunk of
tests of GR itself. It is expected that quantum corrections to GR may add
additional higher powers of $R$ to the Einstein-Hilbert action, or more
generally, modifying the action into a generic class of functions of the Ricci
scalar. Testing the fate of the prized equivalence principle, in such modified
theories of gravity, hence become important in order to obtain a more generic
theory of gravitation, and consequently, of gravitating objects. In this
study, it is shown that a Post-Newtonian (PN) expansion of a class of
$f\left(R\right)$ theories lead to a sequence of solutions to the two-body
problem, which follows the equivalence principle (EP) at the Newtonian order,
and generalizes to the ’effacing principle’ at a higher PN order.

##  I Introduction

The equivalence principle remains to this date a cornerstone in an attempt to
find the most general theory of gravitation possible. Testing the principle of
equivalence has almost become synonymous with testing GR at all possible
length scales, as illustrated in [1]. However, it is entirely possible that
the picture is more subtle than that.

The Equivalence principle (EP), being a principle, is not a theory, and its
consolidation from assumption to fact in the scientific community is based on
empirical (but possibly circumstantial) evidence. Hence it becomes all the
more important to keep on trying to find exceptions to the rule in natural
phenomenon to test the limits of a scientific principle. The physical
consequences of EP are found through the predicted effects of Galilean,
Newtonian, and Einsteinian relativity. Namely, the trajectory of an extended
object in a space-time with other gravitational sources, is independent of the
internal gravitational structure of the extended body, till at a scenario
where the ratio of object extendedness to inter-object distances become
significant. Tidal deformations of individual objects then strongly affect the
gravitational field of the external space-time, and structural details of
individual objects begin to play a significantly important role in the
evolutionary dynamics of the equations of motion of such compact objects.

Various tests of the equivalence principle have been attempted, some of which
date quite far back into the past, like the case of John Philoponus in the 6th
century [2], Galileo’s tests in 1610 [3], and Newton’s pendulum experiments in
1680 [4]. However, the most rigorous tests of the equivalence principle in the
pre-modern times were the torsion pendulum experiments of Eötvös in [5] and
following publications. Modern tests of the EP include [6, 7], whereas strong
field tests of the strong EP was performed in [8, 9]

The most extreme test of the equivalence principle conceived so far is through
the analysis of gravitational wave (GW) data. In GW physics, the test of the
equivalence principle, or a consequence of it, is known as the ’effacing
principle’ in literature [10, 11]. The principle can be stated as follows: the
internal gravitational (strong or weak) details of extended compact objects in
a binary system, neutron stars (NS) or black holes (BH), does not influence
the eventual evolution of the trajectory of the individual objects, till at a
very later stage of the orbital evolution. In other words, the deformation of
a compact object due to the presence of another compact object in its
vicinity, although changing the structure of the first object, does not
influence the gravitational field of the external space-time, and the dynamics
of the center of mass (COM) of the massive compact objects remain unaffected
by said deformations. The effacing principle is broken at scenarios where the
average radii of the objects become comparable to the distances between
objects where tidal deformations have significant effects on the orbital
dynamics. Information about the equations of motion of a binary system, for
example, comes to an observer at asymptotic infinity through the phase
evolution of GWs. Studying such data, one can estimate the efficacy of GR in
terms of obtaining an accurate relationship between orbital dynamics and GW
data.

In terms of testing the principle, if objects move on geodesics of an external
gravitational field, then the rate at which they rotate and fall around
(towards) each other is fixed only by their masses and orbital separation.
Perturbatively solving GR field equations from a Newtonian/Keplerian initial
data leads to a sequence of solutions that obey the effacing principle till
2.5 Post Newtonian order [10], for two non-spinning or slowly spinning compact
objects in a binary system. That is, their equations of motion are dependent
only on their masses and distance between their center of masses, till at a
much later stage of the orbital evolution.

However, consequences of the effacing principle might not be unique to GR. To
show this, we take a class of theories which are generalizations of the
Einstein-Hilbert action of GR, collectively called as $f\left(R\right)$
theories of gravity [12, 13], and expand an auxiliary field about a Minkowski
space-time using the Post-Newtonian formalism. The action $S$ from which the
field equations of $f\left(R\right)$ theories of gravity are obtained is given
by the following 4-integral over a manifold $\mathscr{M}$

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{\mathscr{M}}\,\sqrt{-g}\,\,d^{4}x\,\left[\frac{f\left(R\right)}{2\,\kappa^{2}}\,+\,T_{m}\right]$
(1)

where $f\left(R\right)$ is a generic function of the Ricci scalar $R$, and
$\sqrt{-g}\,d^{4}x$ is the covariant infinitesimal volume element. $T_{m}$ is
the classical matter action. The above reduces to the familiar Einstein-
Hilbert action for $f\left(R\right)\,=\,R$, and the field equations of GR can
hence be recovered. Geometrized units, $c\,=\,G\,=\,1$ will be used in this
study, which implies: $\kappa^{2}\,=\,8\,\pi$. The notations of [14] will be
followed except a few changes in variable and index labeling.

In this article, the two-body problem in $f\left(R\right)$ theories of gravity
will be solved for, using an initial data on the metric and matter that is
asymptotically Newtonian/Keplerian A relaxed form of the field equations will
be used, similar to the Landau-Lifshitz formulation of GR in [15,
*LANDAU1975295, *LANDAU1975259], leading to a post-Newtonian sequence of
solutions, as reviewed in Sec. II. The solutions will be expressed as
functions of multipole moments over two compact sources, and it will be
explicitly shown that the equations of motion of the two sources are dependent
only on the monopole moments of the sources, and the distance between their
’center of masses’.

In Sec. III a method to relax the field equations of $f\left(R\right)$
theories will be illustrated. Using the formalism of Sec. II ”gauge invariant”
modification to the 6 PN equations of motion of GR due to $f\left(R\right)$
will be derived. It will also be shown that even though deviations from
general relativistic equations of motion exist, consequences of the
equivalence (and consequently the effacing) principle appear in the equations
of motion of two compact objects in $f\left(R\right)$ theories of gravity, and
is almost identical to GR.

##  II The Equivalence and the Effacing principle in General Relativity

###  II.1 A coordinate system adapted for the internal problem

An asymptotic series of solutions based on Newtonian initial data can be found
as a function of multipole moments of classical matter using a formalism
developed by Schutz and Futamase [18], by using two co-moving and scaled ’body
zone’ coordinate systems for each of the constituents of the binary, compared
to an observer at Minkowskian asymptotic infinity.

The formalism involves defining a scaled time coordinate $s$ (different from
$\tau$ as was used in [19]), with respect to the time coordinate $t$ used by
an asymptotic observer at rest, with respect to the center of mass of the
binary system,

$\displaystyle s$ $\displaystyle=$ $\displaystyle\epsilon\,t,$ (2)

which is well suited for the internal problem of the two bodies, given they
are non-spinning or slowly spinning, as was found in [18]. $\epsilon$ is the
post Newtonian parameter given by $\frac{\left|v_{orb}\right|}{c}$, with
$\left|v_{orb}\right|$ being the absolute value of the orbital velocity, and
$c$ being the speed of light (unity in the current study). The re-
parameterized time $s$ is also taken as the proper time of a co-moving
observer along either sources.

A transformation into the ’body zone coordinate system’ involves first setting
up of a Fermi-Walker coordinate system along the center of mass of two compact
objects [19, 20, 14], whose centers of masses follow two world lines, denoted
by 3-vectors $z_{L}^{i}\left(s\right)$, for $L\,=\,1,\,2$. Under such a
transformation, various components of a contravariant tensor, as defined
initially in the asymptotic observer’s frame (denoted in primed lowercase
Latin), transform into the co-moving frame as follows

$\displaystyle T_{L}^{s\,s}$ $\displaystyle=$ $\displaystyle T_{L}^{s\,s}$ (3)
$\displaystyle T_{L}^{s\,i}$ $\displaystyle=$
$\displaystyle\epsilon^{2}T_{L}^{s\,i^{\prime}}+v_{L}^{i}T_{L}^{s\,s}$ (4)
$\displaystyle T_{L}^{i\,j}$ $\displaystyle=$
$\displaystyle\epsilon^{4}T_{L}^{i^{\prime}\,j^{\prime}}+2\epsilon^{2}v_{L}^{\left(i\right.}T_{L}^{\left.j^{\prime}\right)s}+v_{L}^{i}\,v_{L}^{j}T_{L}^{s\,s}$
(5)

The ’smallness’ parameter $\epsilon$ is the ratio of the average orbital
velocities and the speed of light (unity in the current study), as agreed upon
by an observer at asymptotic infinity, who is static with respect to the
center of mass of the two bodies. Then the 3-velocity of the body
$v_{L}^{i}\left(\,=\,\frac{dz^{i}_{L}}{ds}\right)$, as measured by an external
asymptotic observer, is of order unity.

The next step involves defining a specially scaled spatial coordinate system
for the internal problem. For compact objects whose mass-radii ratios approach
unity, given the orbital separation is held fixed, both mass and radius scale
as $\epsilon^{2}$, if one were to utilize a Newtonian/Keplerian initial data.
Consequently, mass-energy densities scale as $\epsilon^{-4}$. Hence, the
internal problem can be solved in coordinate systems $X_{L}^{i}$ (as defined
in [19, 20, 14]) that is scaled by the parameter $\epsilon^{2}$, such that

$\displaystyle X_{L}^{i}$ $\displaystyle\equiv$
$\displaystyle\frac{x^{i}-z_{L}^{i}\left(s\right)}{\epsilon^{2}}.$ (6)

Under the above transformation, various contravariant components of the
energy-momentum tensor $T^{\mu\nu}$ scale as follows in the body zone

$\displaystyle T_{L}^{s\,s}$ $\displaystyle=$
$\displaystyle\mathscr{O}\left(\epsilon^{-2}\right)$ (7) $\displaystyle
T_{L}^{s\,i^{\prime}}$ $\displaystyle=$
$\displaystyle\mathscr{O}\left(\epsilon^{-4}\right)$ (8) $\displaystyle
T_{L}^{i^{\prime}\,j^{\prime}}$ $\displaystyle=$
$\displaystyle\mathscr{O}\left(\epsilon^{-8}\right)$ (9)

which completes the set of transformations necessary to solve the internal
problem. It is to be noted that, as seen in the special coordinate systems
(the body zone coordinates), the radii of either compact objects are unity.

###  II.2 Equations of motion

Using the ’special’ coordinate system defined in the previous section, a
sequence of solutions can be found for the internal problem. The relaxed form
of the field equations of GR, with a matter source tensor $T^{\mu\nu}_{m}$, in
the asymptotically Minkowskian observer’s coordinate system are given as
follows

$\displaystyle\partial_{\alpha\,\beta}H^{\mu\alpha\nu\beta}$ $\displaystyle=$
$\displaystyle
16\,\pi\,\left(-g\right)\left(T^{\mu\nu}_{m}+t_{LL}^{\mu\nu}\right)$ (10)
$\displaystyle H^{\mu\alpha\nu\beta}$ $\displaystyle=$
$\displaystyle\mathfrak{g}^{\alpha\beta}\mathfrak{g}^{\mu\nu}-\mathfrak{g}^{\alpha\nu}\mathfrak{g}^{\beta\mu}$
(11)

$\mathfrak{g}^{\mu\nu}$ is the square root of the determinant weighed
contravariant metric $\left(\sqrt{-g}\,g^{\mu\nu}\right)$. In this picture of
GR, the dynamical variable is not the metric $g^{\mu\nu}$, but the
contravariant metric density $\mathfrak{g}^{\mu\nu}$ (also known as the gothic
metric in PN literature), which propagates on a Minkowski background
$\eta^{\mu\nu}$. An apparent separation of the metric of a maximally symmetric
space-time (corresponding to the tensor $H^{\mu\alpha\nu\beta}$), from the
part of the metric that leads to curvature of space-time (the effective
energy-momentum pseudo-tensor of gravitation $t^{\mu\nu}_{LL}$), is possible
in this formalism. The exact form is given by

$\displaystyle
16\,\pi\,\left(-g\right)\,t_{\mathrm{LL}}^{\alpha\beta}=\partial_{\lambda}\mathfrak{g}^{\alpha\beta}\partial_{\mu}\mathfrak{g}^{\lambda\mu}-\partial_{\lambda}\mathfrak{g}^{\alpha\lambda}\partial_{\mu}\mathfrak{g}^{\beta\mu}\,+\,\frac{1}{2}\mathfrak{g}^{\alpha\beta}\mathfrak{g}_{\lambda\mu}\partial_{\rho}\mathfrak{g}^{\lambda\nu}\partial_{\nu}\mathfrak{g}^{\mu\rho}-\mathfrak{g}^{\alpha\lambda}\mathfrak{g}_{\mu\nu}\partial_{\rho}\mathfrak{g}^{\beta\nu}\partial_{\lambda}\mathfrak{g}^{\mu\rho}$
$\displaystyle\,-\,\mathfrak{g}^{\beta\lambda}\mathfrak{g}_{\mu\nu}\partial_{\rho}\mathfrak{g}^{\alpha\nu}\partial_{\lambda}\mathfrak{g}^{\mu\rho}\,+\,\mathfrak{g}_{\lambda\mu}\mathfrak{g}^{\nu\rho}\partial_{\nu}\mathfrak{g}^{\alpha\lambda}\partial_{\rho}\mathfrak{g}^{\beta\mu}\,+\,\frac{1}{8}\left(2\mathfrak{g}^{\alpha\lambda}\mathfrak{g}^{\beta\mu}\,-\,\mathfrak{g}^{\alpha\beta}\mathfrak{g}^{\lambda\mu}\right)\left(2\mathfrak{g}_{\nu\rho}\mathfrak{g}_{\sigma\tau}\,-\,\mathfrak{g}_{\rho\sigma}\mathfrak{g}_{\nu\tau}\right)\partial_{\lambda}\mathfrak{g}^{\nu\tau}\partial_{\mu}\mathfrak{g}^{\rho\sigma}$
(12) $\displaystyle\mathfrak{g}_{\mu\nu}\equiv\frac{g_{\mu\nu}}{\sqrt{-g}}$
(13)

which has at most one derivative of $\mathfrak{g}^{\alpha\beta}$. The purpose
of the metric $g$ in this formalism is solely to raise and lower indices of
the auxiliary variable $\mathfrak{g}^{\alpha\beta}$. It is to be noted that
the above definition is neither covariant nor unique. A covariant proposition
for the Landau-Lifshitz energy-momentum pseudo-tensor can be found in [21].

From the old days of Schwarzschild [22], the relaxed form of the field
equations, or the Landau-Lifshitz formalism, helps one to find exact solutions
of the highly non-linear field equations of GR, as was also explained clearly
in [23].

With regards to observable effects, the biggest accomplishment of the relaxed
formalism has been to obtain a series of potentials known as the Post-
Newtonian expansion of GR in literature, as in [24, 25], and later
systematized in [26, *Chandra1967, *Chandra1969I, *Chandra1969II,
*Chandra1970]. It involves defining tensor potentials $\mathfrak{h}^{\mu\nu}$,
given by

$\displaystyle\mathfrak{h}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\eta^{\mu\nu}-\mathfrak{g}^{\mu\nu}$ (14)

along with a constraint

$\displaystyle\partial_{\mu}\mathfrak{h}^{\mu\nu}$ $\displaystyle=$
$\displaystyle 0,$ (15)

which restricts the possible class of coordinates (or gauges), in which one
attempts to find solutions for $\mathfrak{h}^{\mu\nu}$. In PN and GW
literature, the constraint is known as the harmonic gauge, in which the
coordinate chart $x^{\mu}$ follow four massless and homogeneous wave equations

$\displaystyle\Box\,x^{\mu}$ $\displaystyle=$ $\displaystyle 0$ (16)

where $\Box$ is the d’Alembert operator or
$\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$. Such a class of coordinates are
consequently known as harmonic coordinate systems, and the dynamical equations
for $\mathfrak{h}^{\mu\nu}$ take the following form

$\displaystyle\Box\mathfrak{h}^{\mu\nu}$ $\displaystyle=$
$\displaystyle-16\,\pi\,\Lambda^{\mu\nu}$ (17) $\displaystyle\Lambda^{\mu\nu}$
$\displaystyle=$
$\displaystyle\Theta^{\mu\nu}+\partial_{\alpha\,\beta}\chi^{\mu\nu\alpha\beta}$
(18) $\displaystyle\Theta^{\mu\nu}$ $\displaystyle=$
$\displaystyle\left(-g\right)\left(T^{\mu\nu}\,+\,t_{LL}^{\mu\nu}\,+\,t_{H}^{\mu\nu}\right)$
(19) $\displaystyle\chi^{\mu\nu\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{16\,\pi}\left(\mathfrak{h}^{\alpha\nu}\mathfrak{h}^{\beta\mu}-\mathfrak{h}^{\alpha\beta}\mathfrak{h}^{\mu\nu}\right)$
(20)

Eq. (17) can be solved for $\mathfrak{h}^{\mu\nu}$ in a series of multipole
moments of the energy-momentum pseudo-tensor $\Lambda^{\mu\nu}$. It is to be
noted that the RHS of Eq. (20) manifests as a total derivative term in the
field equations and can be transformed into boundary terms using Gauss’s law.
Choice of the boundary condition fixes the form of the Green’s function. For
the case in study, the no incoming radiation from past null Minkowskian
infinity boundary condition will be used, which was given in Sec. 4.1 of [14]
as follows

$\displaystyle\lim\limits_{s\,\rightarrow\,r\,,\,r\,\rightarrow\,\infty}\,\,\left[\frac{\partial}{\partial\,r}\left(r\,\mathfrak{h}^{\mu\nu}\right)\,+\,\frac{\partial}{\partial\,s}\left(r\,\mathfrak{h}^{\mu\nu}\right)\right]$
$\displaystyle=$ $\displaystyle 0$ (21)

where $r$ is the radial coordinate distance from the COM of a compact source
to a field point. Enforcing the condition (21), the effect of the terms in the
RHS of Eq. (20) on the solution sequence $\mathfrak{h}^{\mu\nu}$ vanishes.
Various components of the tensor potentials $\mathfrak{h}^{\mu\nu}$ till
$\mathscr{O}\left(\epsilon^{6}\right)$, in the body zone coordinate system,
are given as follows [19, 20, 14]

$\displaystyle\mathfrak{h}_{\mathscr{B}}^{s\,s}$ $\displaystyle=$
$\displaystyle
4\,\epsilon^{4}\,\sum_{L=1,\,2}\,\left(\frac{P^{s}_{L}}{\left|Z^{i}_{L}\right|}+\epsilon^{2}\,\frac{D^{k}_{L}\,Z^{k}_{L}}{\left|Z^{i}_{L}\right|^{3}}\right)+\mathscr{O}\left(\epsilon^{8}\right)$
(22) $\displaystyle\mathfrak{h}_{\mathscr{B}}^{s\,i}$ $\displaystyle=$
$\displaystyle
4\,\epsilon^{4}\sum_{L=1,\,2}\,\left(\frac{P_{L}^{i}}{\left|Z^{i}_{L}\right|}+\epsilon^{2}\,\frac{J_{L}^{k\,i}Z_{L}^{k}}{\left|Z^{i}_{L}\right|^{3}}\right)+\mathscr{O}\left(\epsilon^{8}\right)$
(23) $\displaystyle\mathfrak{h}_{\mathscr{B}}^{i\,j}$ $\displaystyle=$
$\displaystyle
4\,\epsilon^{2}\,\sum_{L=1,\,2}\left(\frac{Y_{L}^{i\,j}}{\left|Z^{i}_{L}\right|}+\epsilon^{2}\,\frac{Y_{L}^{k\,i\,j}\,Z_{L}^{k}}{\left|Z^{i}_{L}\right|^{3}}\right)+\mathscr{O}\left(\epsilon^{6}\right)$
(24)

where $Z^{i}_{L}\left(s\right)\,=\,x^{i}\,-\,z^{i}_{L}\left(s\right)$ is the
displacement 3-vector of the field point to the COM of the source. In the
current formalism, the above components of the tensor potential are sufficient
to obtain equations of motion of the two sources till the first PN order. The
monopole and dipole moments, respectively, of various components of
$\Lambda^{\mu\nu}$ are defined as follows

$\displaystyle P_{L}^{s}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,\int_{\mathscr{B}_{L}}d^{3}X_{L}\,\Lambda^{s\,s}$
(25) $\displaystyle P^{i}_{L}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,\int_{\mathscr{B}_{L}}d^{3}X_{L}\,\Lambda^{s\,i}$
(26) $\displaystyle Y^{i\,j}_{L}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,\int_{\mathscr{B}_{L}}d^{3}X_{L}\,\Lambda^{i\,j}$
(27) $\displaystyle D^{i}_{L}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,\int_{\mathscr{B}_{L}}d^{3}X_{L}\,\Lambda^{s\,s}\,X^{i}_{L}$
(28) $\displaystyle J^{i\,j}_{L}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,\int_{\mathscr{B}_{L}}d^{3}X_{L}\,\Lambda^{s\,i}X^{j}_{L}$
(29) $\displaystyle Y^{i\,j\,k}_{L}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,\int_{\mathscr{B}_{L}}d^{3}X_{L}\,\Lambda^{i\,j}\,X^{k}_{L}$
(30)

where the various pre-factors of $\epsilon^{n}$ for $n\,=\,2,\,4,\,8$ for the
various multipole moment definitions have been removed, compared to [19, 20,
14], in order to avoid any confusion regarding the scaling of various
components of a rank two tensor $\Lambda^{s\,s}$ under the transformations (2)
and (6).

The equations of motion at various PN orders come from the definition of a
quasi-local four-momenta of the COM of each of the bodies in the binary as
defined in [19, 20, 14, 31] as

$\displaystyle P^{\mu}_{L}\left(s\right)$ $\displaystyle=$
$\displaystyle\epsilon^{2}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\Lambda^{s\,\mu},$
(31)

and the conservation law

$\displaystyle\Lambda^{\mu\nu}_{\quad\,,\,\nu}$ $\displaystyle=$
$\displaystyle 0,$ (32)

from which one obtains the evolution equation of the four-momenta
$P^{\mu}_{L}$ as surface integrals of the energy-momentum pseudo-tensor, over
the boundary of the body zone of the $L^{th}$ object, as was obtained in [19,
20, 14]:

$\displaystyle\frac{dP^{\mu}_{L}}{ds}$ $\displaystyle=$
$\displaystyle-\epsilon^{-4}\,\oint_{\partial\mathscr{B}_{L}}\,dS_{k}\,\Lambda^{k\,\mu}+\epsilon^{-4}\,v^{k}_{L}\,\oint_{\partial\mathscr{B}_{L}}\,dS_{k}\,\Lambda^{s\,\mu},$

where $dS_{k}$ is an infinitesimal unit vector normal to the 2-sphere
$\partial\mathscr{B}_{L}$ that is the boundary of the body zone
$\mathscr{B}_{L}$. The 3-momentum vs 3-velocity relationship is given by

$\displaystyle P^{i}_{L}$ $\displaystyle=$ $\displaystyle
P^{s}_{L}\,v^{i}_{L}+Q^{i}_{L}+\mathscr{O}\left(\epsilon^{2}\right)$ (34)
$\displaystyle Q^{i}_{L}$ $\displaystyle=$
$\displaystyle\epsilon^{-4}\oint_{\partial\mathscr{B}_{L}}dS_{k}\left(\Lambda^{s\,k}-v_{L}^{k}\,\Lambda^{s\,s}\right)X_{L}^{i}$
(35)

Newtonian like equations of motion involve evolution equations of
3-velocities, which are the 3-accelerations. They were given in [19, 20, 14]
as

$\displaystyle P^{s}_{L}\,\frac{dv^{i}_{L}}{ds}$ $\displaystyle=$
$\displaystyle-\epsilon^{-4}\oint_{\partial\mathscr{B}_{L}}dS_{k}\,\Lambda^{k\,i}+\epsilon^{-4}\,v^{k}_{L}\oint_{\partial\mathscr{B}_{L}}dS_{k}\,\Lambda^{s\,i}$
(36)
$\displaystyle+\epsilon^{-4}\,v^{i}_{L}\left(\oint_{\partial\mathscr{B}_{L}}dS_{k}\,\Lambda^{k\,s}-v_{L}^{k}\oint_{\partial\mathscr{B}_{L}}dS_{k}\,\Lambda^{s\,s}\right)$
$\displaystyle-\frac{dQ^{i}_{L}}{ds},$

If one chooses a definition of the mass $M_{L}$ of the $L^{th}$ body as

$\displaystyle M_{L}$ $\displaystyle=$
$\displaystyle\lim\limits_{\epsilon\,\rightarrow\,0}\,P^{s}_{L},$ (37)

then the equations of motion of the first body, for example, till first PN
order were obtained in [19, 20, 14] as

$\displaystyle M_{1}\,\frac{dv_{1}^{i}}{ds}$ $\displaystyle=$
$\displaystyle-\frac{M_{1}\,M_{2}}{r_{1\,2}^{\,2}}\,n^{i}\,+\,\epsilon^{2}\,\frac{M_{1}\,M_{2}}{r_{1\,2}^{2}}\,\left[\left(-v_{1}^{2}-2\,v_{2}^{2}\right.\right.$
$\displaystyle\left.\left.+\frac{3}{2}\,\left(\hat{n}\,\cdot\textbf{v}_{2}\right)^{2}+4\,\left(\textbf{v}_{1}\,\cdot\,\textbf{v}_{2}\right)+\frac{5\,M_{1}}{r_{1\,2}}+\frac{4\,M_{2}}{r_{1\,2}}\right)\,n^{i}\right.$
$\displaystyle\left.\,+\,\left\\{4\,\left(\hat{n}\,\cdot\,\textbf{v}_{1}\right)\,-3\,\left(\hat{n}\,\cdot\,\textbf{v}_{2}\right)\right\\}\,\left(v_{1}^{i}-v_{2}^{i}\right)\right]\,+\,\mathscr{O}\left(\epsilon^{4}\right)$

where $n^{i}\,\left(\,\equiv\,\hat{n}\right)$ is a unit vector pointing from
the COM of the first body zone to the COM of the second body zone, and
$r_{1\,2}$ is the distance between the COM of the two body zones or the
orbital separation. Mass multipole moments of $\Lambda^{\mu\nu}$, as seen by
an observer at rest (w.r.t. the COM of the binary system) at asymptotic
infinity, go as $\left(\text{mass}\right)\,\times\,\left(\text{average radius
of mass distribution}\right)^{\ell}$. Using Newtonian initial data, mass
multipole objects appearing in the equations of motion of the binary (36)
scale as $\mathscr{O}\left(\epsilon^{2\,\ell\,+\,2}\right)$. In terms of the
classical notion of force on a particle,

$\displaystyle F^{i}$ $\displaystyle=$ $\displaystyle
F^{i}_{N}\,+\,\epsilon^{2}\,F^{i}_{1PN}\,+\,\mathscr{O}\left(\epsilon^{4}\right)$
(39)

corresponding to the first and second term of the RHS of Eq. (II.2),
respectively. If one simply uses the first term of (II.2), as a crude
approximation for the trajectory of two slow spinning compact objects in the
early inspiral phase, one obtains the Newtonian force on the first body. The
dynamics of the first body then is independent of its own mass and depends
only on the Newtonian potential generated by the second body, which is a
consequence of the equivalence principle. The generalization of the
consequences of the equivalence principle till first PN order comes through
the effacing principle, where the internal/structural details of either body,
and their effect on the external gravitational field, do not factor into their
trajectories around/towards each other. The compact objects still behave like
point particles moving along the geodesics of an external gravitational field.
Self force effects through the appearance of velocity dependent terms appear
at the first PN order. Although from the first PN order onward the dynamics of
the first body zone is dependent on the mass $M_{1}$ enclosed in it, as was
defined in Eq. (37), there exist no multipole objects other than the monopole
moments of $\Lambda^{\mu\nu}$, essentially making the two compact objects
behave like massive but point particles (not to be confused with test masses).

The effects of spin in the equations of motion are dependent on the scaling of
the current multipole moments of $\Lambda^{\mu\nu}$. In the slow rotation
approximation, the internal velocities scale as
$\mathscr{O}\left(\epsilon\right)$. Hence current multipole moments go as
$\left(\text{mass}\right)\,\times\,\left(\text{average radius of mass
distribution}\right)^{\ell}\,\times\,\left(\text{velocity of internal
motion}\right)$, which is $\mathscr{O}\left(\epsilon^{2\,\ell\,+\,3}\right)$,
implying that the spin-orbit coupling force is of the form
$\left(\text{mass}\right)\,\times\,\left(\text{orbital
velocity}\right)\,\times\,\left(\text{spin}\right)$, or
$\mathscr{O}\left(\epsilon^{2\,\ell\,+\,4}\right)$. It is also to be noted
that the scalings change when the velocity of internal motion cannot be
ignored for rapidly rotating constituents of the binary. The time scaling in
Eq. (2), then changes to $s\,\equiv\,\epsilon^{-2}\,t$, as was described in
[19].

The above described consequence of GR is a way one may choose to interpret the
effects of the equivalence assumption. But this kind of a consequence is not
unique to GR, as will be shown in the following section.

##  III First PN equations of motion in $f(R)$ gravity and the Effacing
principle

###  III.1 Relaxing the $f\left(R\right)$ field equations

The field equations of $f\left(R\right)$ theories of gravity can be found by
varying the action in Eq. (1). It is however, for the sake of understanding
the physical consequences of such theories, better to frame the field
equations in an Einsteinian way; such that possible observable deviations from
GR can be expressed as an effective energy-momentum tensor, that is different
from the ordinary matter energy-momentum tensor. In the latter way, the
contravariant field equations are expressed as [32]

$\displaystyle
G^{\mu\nu}=\frac{8\,\pi}{f^{\prime}}\,\left(T_{m}^{\mu\nu}\,+\,T^{\mu\nu}_{eff}\right)$
(41) $\displaystyle
3\,g^{\mu\nu}\,\nabla_{\mu}\,\nabla_{\nu}\,f^{\prime}\,+\,f^{\prime}\,R\,-\,2\,f=8\,\pi\,T_{m}$
$\displaystyle T^{\mu\nu}_{eff}$ $\displaystyle\equiv$
$\displaystyle\nabla^{\mu}\,\nabla^{\nu}\,f^{\prime}\,+\,\frac{g^{\mu\nu}}{2}\,\left(f\,-\,R\,f^{\prime}\right)\,-\,g^{\mu\nu}\,g^{\alpha\beta}\,\nabla_{\alpha}\,\nabla_{\beta}\,f^{\prime}$
$\displaystyle f^{\prime}\left(R\right)$ $\displaystyle\equiv$ $\displaystyle
f^{\prime}\,=\,\frac{df\left(R\right)}{dR}\qquad\,T_{m}\,=\,g_{\alpha\beta}\,T_{m}^{\alpha\beta}\qquad\nabla^{\mu}\equiv
g^{\mu\nu}\,\nabla_{\nu}$

where the subscript $m$ has been utilized to distinguish the classical matter
energy-momentum tensor from the effective one. $G^{\mu\nu}$ is the Einstein
tensor of GR, and $\nabla_{\mu}$ is the covariant derivative associated with
the metric $g^{\mu\nu}$.

The above system of partial differential equations are extremely non-linear in
the metric and the scalar field $f^{\prime}$, and needs some form of relaxing
before a solution, similar to the PN expansion in GR, can be generated as
functions of multipole moments of an effective source given by the RHS of Eq.
(41).

There has been a number of analysis on $f\left(R\right)$ theories of gravity.
In the cosmological context and local gravity constraints, a review can be
found in [33, 34, 35]. Early solutions for spherically symmetric distributions
of fluids can be found in [36, 37, 38, 39, 40]. An analysis using null tetrads
to obtain conditions for deviations to a curvature invariant object (Weyl
scalar) $\Psi_{2}$ was done in [41]. Perturbations about a completely
homogeneous and isotropic space-time was performed in [42], whereas
perturbation studies about a spherically symmetric space-time using the
standard metric PN expansion (without using the gothic metric density) was
done in [43]. Stability of solutions about perturbations in general space-
times for various deviation parameters was shown in [44].

In the spirit of the Landau-Lifshitz formalism, an auxiliary metric
$\mathfrak{\tilde{g}}^{\mu\nu}$, similar to the gothic metric
$\mathfrak{g}^{\mu\nu}$ of GR, but different from [42, 43] can be defined as

$\displaystyle\mathfrak{\tilde{g}}^{\mu\nu}$ $\displaystyle=$ $\displaystyle
f^{\prime}\left(R\right)\,\sqrt{-g}\,g^{\mu\nu}$ (44)

which will be the dynamical variable in the current study, propagating in a
Minkowski background $\eta^{\mu\nu}$. It is to be noted that, under
linearization of the field equations of $f\left(R\right)$ theories of gravity
with a well defined polynomial expansion about $R\,=\,0$, the RHS of Eq. (44)
reduces to the transverse-traceless metric deviation tensor $\bar{h}^{\mu\nu}$
in [40].

Using the new definition, the field equations (41) can be written as follows

$\displaystyle\partial_{\alpha\,\beta}H^{\alpha\mu\beta\nu}$ $\displaystyle=$
$\displaystyle
16\,\pi\,\left(-g\right)\,f^{\prime}\left(R\right)\,\left(T_{m}^{\mu\nu}\,+\,t_{eff}^{\mu\nu}\,+\,t^{\mu\nu}_{LL}\right),$
(45) $\displaystyle 16\,\pi\,\left(-g\right)\,f^{\prime}\,t^{\mu\nu}_{eff}$
$\displaystyle=$
$\displaystyle\sqrt{-g}\,\left(f\,-\,R\,f^{\prime}\right)\,\tilde{\mathfrak{g}}^{\mu\nu}\,+\,\frac{3}{\left(f^{\prime}\right)^{2}}\,\left(\tilde{\mathfrak{g}}^{\mu\alpha}\,\tilde{\mathfrak{g}}^{\nu\beta}\,-\,\frac{1}{2}\,\tilde{\mathfrak{g}}^{\mu\nu}\,\tilde{\mathfrak{g}}^{\alpha\beta}\right)\,\partial_{\alpha}f^{\prime}\,\partial_{\beta}f^{\prime}$
(46) $\displaystyle f$ $\displaystyle\equiv$ $\displaystyle
f\left(R\right)\,;\,\qquad\qquad f^{\prime}\,\equiv\,f^{\prime}\left(R\right)$
(47)

where the forms of $H^{\alpha\mu\beta\nu}$ and $t^{\mu\nu}_{LL}$ remain the
same as in GR (with $\mathfrak{g}$ replaced by $\tilde{\mathfrak{g}}$) as was
defined in Eqs. (11) and (II.2), whereas $t^{\mu\nu}_{eff}$ is seen to be a
quadratic function of only first derivatives of $f^{\prime}\left(R\right)$.
The functional forms of $f\left(R\right)$ and $f^{\prime}\left(R\right)$ can
be written as a power series about $R\,=\,0$ as

$\displaystyle f\left(R\right)$ $\displaystyle=$ $\displaystyle
R\,+\,\frac{f^{\prime\prime}\left(0\right)}{2}\,R^{2}$ (48) $\displaystyle
f^{\prime}\left(R\right)$ $\displaystyle=$ $\displaystyle
1\,+\,f^{\prime\prime}\left(0\right)\,R$ (49)

where the coefficients of expansion $f^{\prime}\left(0\right)$ is taken to be
unity to recover GR at the $R\,=\,0$ limit, and
$f^{\prime\prime}\left(0\right)$ will be taken to be negative for the course
of this article, following [40, 36, 37, 39]. The negative choice for the
constant $f^{\prime\prime}\left(0\right)$ also frames Eq. (41), after making
the choices (48) and (49), into a Klein-Gordon equation that is satisfied by
the massive scalar field.

Perturbing the configuration $\left(\mathfrak{\tilde{g}},\,R\right)$ about a
Minkowski space-time $\eta^{\mu\nu}$, Eq. (45) can now be perturbatively
solved as functions of multipole moments of net energy-momentum pseudo-tensors
$\left(\tilde{\Lambda}^{\mu\nu},\,\tilde{\Lambda}\right)$. One can define a
tensor potential $\tilde{\mathfrak{h}}^{\mu\nu}$, similar to
$\mathfrak{h}^{\mu\nu}$ of GR, that propagates on a Minkowski space-time
$\eta^{\mu\nu}$

$\displaystyle\mathfrak{\tilde{h}}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\eta^{\mu\nu}\,-\,\mathfrak{\tilde{g}}^{\mu\nu},$ (50)

and use a conformal-Lorenz gauge condition (referred to as the ”Lorentz” gauge
condition in [45], and as conformal-harmonic condition in [46])

$\displaystyle\mathfrak{\tilde{h}}^{\mu\nu}_{\quad\,,\,\nu}$ $\displaystyle=$
$\displaystyle 0.$ (51)

Since there is the Ricci scalar manifesting as a scalar field, one needs to
define a dynamical variable that propagates on a Ricci flat background, and
has a ”proper” scaling relationship with the trace of the classical energy-
momentum tensor (which may comprise of a pair of slowly spinning fluids,
gravitational mass monopoles, or a combination of both). Noticing that
$\Theta^{\mu\nu}$ has an overall factor of $\left(-g\right)$ in the RHS of Eq.
(19), such that $\Box\,\mathfrak{h}^{\mu\nu}$ is related to the classical
energy-momentum tensor $T_{m}^{\mu\nu}$ with an overall $\left(-g\right)$
factor. One may choose to structure the trace of the field equations (41) in a
manner that connects the d’Alembert operated scalar density, say
$\Box\,\mathfrak{R}$, to the trace of $T_{m}^{\mu\nu}$ with an overall
$\left(-g\right)$ factor.

In the action, the square root of the negative of metric determinant
multiplied with the Ricci scalar is akin to the energy-density (if one uses
geometrized units). Hence the scalar variable, like the tensor variable
$\mathfrak{g}^{\mu\nu}$ in GR and $\tilde{\mathfrak{g}}^{\mu\nu}$ in
$f\left(R\right)$, is changed to be the following

$\displaystyle\mathfrak{R}$ $\displaystyle=$ $\displaystyle\sqrt{-g}\,R$ (52)

One obtains the relaxed dynamics of $\mathfrak{\tilde{h}}^{\mu\nu}$ as

$\displaystyle\Box\mathfrak{\tilde{h}}^{\mu\nu}$ $\displaystyle=$
$\displaystyle 16\,\pi\,\tilde{\Lambda}^{\mu\nu}$ (53)
$\displaystyle\tilde{\Lambda}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\left(-g\right)\,\left[T^{\mu\nu}_{m}\,\left(1\,-\,\frac{f^{\prime\prime}\,\,\mathfrak{R}}{\sqrt{-g}}\right)\,+\,t_{LL}^{\mu\nu}\,+\,t_{H}^{\mu\nu}\,+\,t_{eff}^{\mu\nu}\right].$
(54)

with the following conservation law being satisfied by
$\tilde{\Lambda}^{\mu\nu}$

$\displaystyle\tilde{\Lambda}^{\mu\nu}_{\,\,\,,\,\mu}$ $\displaystyle=$
$\displaystyle 0\,\,\,.$ (55)

$t_{eff}^{\mu\nu}$ is comprised of various products of
$\tilde{\mathfrak{h}}^{\mu\nu}$, $\mathfrak{R}$, and first derivatives of
$\mathfrak{\tilde{h}}^{\mu\nu}$ and $\mathfrak{R}$; whose truncated form till
the quadratic order of $f^{\prime\prime}\left(0\right)$ is given as follows

$\displaystyle
16\,\pi\,\left(-g\right)\,t^{\mu\nu}_{eff}=\,\frac{f^{\prime\prime}\,\,\mathfrak{R}^{2}}{\sqrt{-g}}\,\,\eta^{\mu\nu}\,+\,\frac{3\,\left(f^{\prime\prime}\right)^{2}}{4\,\left(-g\right)}\,\left[2\,\left(\eta^{\mu\alpha}\,\eta^{\nu\beta}\,+\,\eta^{\mu\beta}\,\eta^{\nu\alpha}\right.\right.$
$\displaystyle\left.\left.\,-\,\eta^{\mu\nu}\,\eta^{\alpha\beta}\right)\,\mathfrak{R}\,\partial_{\alpha}\mathfrak{R}\,\partial_{\beta}\tilde{\mathfrak{h}}\,+\,\left(\eta^{\mu\alpha}\,\eta^{\nu\beta}\,-\,\frac{1}{2}\,\eta^{\mu\nu}\,\eta^{\alpha\beta}\right)\,\left(4\,\partial_{\alpha}\mathfrak{R}\,\partial_{\beta}\mathfrak{R}\,+\,\mathfrak{R}^{2}\,\partial_{\alpha}\tilde{\mathfrak{h}}\,\partial_{\beta}\tilde{\mathfrak{h}}\right)\right]$
$\displaystyle\,+\,\mathscr{O}\left[\left(\frac{f^{\prime\prime}}{\sqrt{-g}}\right)^{3}\right]$
(56)
$\displaystyle\tilde{\mathfrak{h}}\,\equiv\,\eta_{\mu\nu}\,\tilde{\mathfrak{h}}^{\mu\nu}$
(57)

It is to be noted that under the transformation (2), the covariant Minkowski
metric, using the $\left(+\,,\,-\,,\,-\,,\,-\right)$ signature becomes

$\displaystyle\eta_{\mu\nu}$ $\displaystyle\equiv$ $\displaystyle
diag\left(\epsilon^{-2},\,-1,\,-1,\,-1\right)$ (58)

therefore
$\left|\det\left(\eta_{\mu\nu}\right)\right|\,\equiv\,\left(-g\right)$ scale
as $\epsilon^{-2}$. It is to be noted that the trace of the new (and the old)
gothic metric density deviation $\tilde{\mathfrak{h}}$ scale as $\epsilon^{2}$
under the transformation (2), compared to the individual components of
$\tilde{\mathfrak{h}}^{\mu\nu}$ at the spinless or slowly spinning limit,
which scale as $\epsilon^{4}$. It is important to note the coefficients of
$\left(-g\right)^{-1}$ in the second term of Eq. (56) (inside the square
brackets) as a whole are sub-leading (in powers of $\epsilon$ in the PN
expansion), compared to the first term of Eq. (56) which is in turn sub-
leading with respect to $T_{m}^{\mu\nu}$. Notice that the first appearance of
a non GR term in $\tilde{\Lambda}^{\mu\nu}$ in the $\epsilon$ sequence is a
non-minimal coupling between the Ricci scalar and the classical matter energy-
momentum tensor along with a term quadratic in $\mathfrak{R}$.

$t_{H}^{\mu\nu}$ appears in the PN formulation of GR as well, whose form was
given in Eqs. (18) and (20) through the term
$\chi^{\mu\nu\alpha\beta}_{\quad\,,\,\alpha\beta}$, with
$\mathfrak{h}^{\mu\nu}$ replaced by $\tilde{\mathfrak{h}}^{\mu\nu}$. The total
derivative term allows one to put $t^{\mu\nu}_{H}$ to zero at the boundaries
of integration because of the specific boundary condition choice of no
incoming tensor radiation from past null Minkowskian infinity.

Under the transformations in Eqs. (44) and (52), the trace of the field
equations, that is Eq. (41), reduces to the following

$\displaystyle\Box\,\mathfrak{R}\,+\,\sqrt{-g}\,\gamma^{2}\,\mathfrak{R}=-8\,\pi\,\left(-g\right)\,\gamma^{2}\,\tilde{\Lambda}$
(59)
$\displaystyle\tilde{\Lambda}=\,T_{m}\,+\,f^{\prime\prime}\left(0\right)\,\left\\{\frac{T_{m}\,\mathfrak{R}}{3\,\sqrt{-g}}+\,\frac{1}{8\,\pi\,\left(-g\right)}\left(\frac{\mathfrak{R}^{2}}{3}\,-\,\partial_{\mu}\mathfrak{R}\,\partial^{\mu}\tilde{\mathfrak{h}}\,-\,\frac{1}{2}\,\mathfrak{R}\,\partial_{\mu}\tilde{\mathfrak{h}}^{\alpha\beta}\,\partial^{\mu}\tilde{\mathfrak{h}}_{\alpha\beta}\right.\right.$
$\displaystyle\left.\left.\,-\,\frac{1}{4}\,\mathfrak{R}\,\partial_{\mu}\tilde{\mathfrak{h}}\,\partial^{\mu}\tilde{\mathfrak{h}}\,-\,\frac{1}{2}\,\mathfrak{R}\,\Box\tilde{\mathfrak{h}}\right)\right\\}\,+\,\mathscr{O}\left(\left[f^{\prime\prime}\left(0\right)\right]^{2}\right)$
(60)
$\displaystyle\gamma^{2}\equiv-\frac{1}{3\,f^{\prime\prime}\left(0\right)}$
(61) $\displaystyle\,T_{m}\,=\,\eta_{\mu\nu}\,T_{m}^{\mu\nu}$ (62)

Derivation of the above differential equation from the trace of
$f\left(R\right)$ field equations have been given in Appendix A. The source
side of Eq. (59) has an energy-momentum scalar $\tilde{\Lambda}$ that is
multiplied by a factor of $\left(-g\right)$, similar to the energy-momentum
tensor $\tilde{\Lambda}^{\mu\nu}$ of Eq. (53). It is to be noted that under
the transformation (2), the presence of the factor of $\sqrt{-g}$ puts an
$\epsilon^{-1}$ dependence on the inverse length squared factor $\gamma^{2}$.
This implies that for $\epsilon\,\rightarrow\,0$, the inverse length is
essentially infinite, and the scalar field $\mathfrak{R}$ cannot be excited.
This also implies that the scalar field excitation is more for values of
$\epsilon$ close to unity, that is strong curvature along with high velocities
in the strong field point particle formalism ([14]). This is known as the
Chameleon mechanism or screening effect in the literature of $f\left(R\right)$
and scalar-tensor theories.

At this point one can impose the boundary conditions that have no radiative
degrees of freedom (time dependent), tensor or scalar, coming from past null
Minkowskian infinity or from any other parts of the space-time. This
assumption implies that the objects in question are isolated objects at
$s\,=\,0$, unaffected by other events in the universe. While this assumption
gets rid of tensor/scalar radiation, one still has to consider non-radiative
and time independent solutions of the homogeneous Klein-Gordon equation.

###  III.2 Homogeneous solution of the Klein-Gordon equation

The time independent and homogeneous Klein-Gordon equation, being a second
order differential equation of only spatial coordinates, is given by

$\displaystyle\nabla\,\mathfrak{R}_{hom}\,-\,\frac{\gamma^{2}}{\epsilon}\,\,\mathfrak{R}_{hom}\,$
$\displaystyle=$ $\displaystyle\,0$ (63)

where $\nabla$ is the Laplace operator in Minkowski space-time, and the square
root of the negative metric determinant was written in terms of the PN
parameter, as seen after Eq. (58). It is seen that the second term in the LHS
of Eq. (63) is of leading order in an $\epsilon\to 0$ series, and the system
as a whole corresponds to the Klein-Gordon equation of a scalar field with
infinite mass for $\epsilon\to 0$. Multiplying both sides of Eq. (63) with
$\epsilon$, and taking the limit $\epsilon\to 0$, one obtains

$\displaystyle\mathfrak{R}_{hom}$ $\displaystyle=$ $\displaystyle 0$ (64)

implying that the homogeneous solution cannot be excited at the leading order,
or when the separation between a binary is large enough / their relative
velocities are vanishing. However, at larger values of $\epsilon$, the first
term in the LHS of Eq. (63) cannot be ignored, that is, Eq. (64) does not hold
true. In spherical symmetry, for example, Eq. (63) can be written as

$\displaystyle\frac{1}{r^{2}}\,\left[r^{2}\,\left(\mathfrak{R}_{hom}\right)_{,r}\right]_{,r}\,-\,\frac{\gamma^{2}}{\epsilon}\,\mathfrak{R}_{hom}$
$\displaystyle=$ $\displaystyle 0$ (65)

which has a solution which is both regular at $r\,=\,0$ and $r\,=\,\infty$,
and is given by

$\displaystyle\mathfrak{R}_{hom}$ $\displaystyle=$
$\displaystyle\frac{C\,\,e^{-\frac{\gamma\,r}{\sqrt{\epsilon}}}}{r}$ (66)

where $C$ is a real constant of integration yet to be fixed. To fix such an
arbitrary integration constant, it is necessary to look for the effect of the
homogeneous solution on the solution of redefined goth metric perturbation, or
$\tilde{\mathfrak{h}}^{\mu\nu}$. Specifically, one needs to find the order of
$\tilde{\mathfrak{h}}^{\mu\nu}$ at which the effect of the scalar homogeneous
solution first appear as a part of the effective source tensor
$t_{eff}^{\mu\nu}$.

###  III.3 The metric deviation in the presence of the homogeneous Ricci
scalar

Since the leading order Ricci deviation is $\mathscr{O}\left(1\right)$ due to
the presence of a homogeneous solution $\mathfrak{R}_{hom}$, it needs to be
included in $\tilde{\Lambda}^{\mu\nu}$ which sources
$\tilde{\mathfrak{h}}^{\mu\nu}$. The solution for
$\tilde{\mathfrak{h}}^{\mu\nu}$ within the boundary condition (21) is given by
the following integral

$\displaystyle\tilde{\mathfrak{h}}^{\mu\nu}\left(s,\,x^{i}\right)$
$\displaystyle=$ $\displaystyle
4\,\left(-g\right)\,\int_{C\left(s,\,x^{k};\,\epsilon\right)}\,d^{3}y\,\frac{\tilde{\Lambda}^{\mu\nu}\left(s\,-\,\epsilon\,\left|x^{i}\,-\,y^{i}\right|,\,y^{i}\right)}{\left|x^{i}\,-\,y^{i}\right|}$
(67)

where the integral has its support in the past light cone of the event
$\left(s,\,x^{i}\right)$, as fixed by the boundary condition, and denoted by
$C\left(s,\,x^{k};\,\epsilon\right)$. Concentrating only on the time-time
component of the above, for example, under the temporal transformation (2),
$\tilde{\mathfrak{h}}^{t\,t}\,=\,\epsilon^{-2}\,\tilde{\mathfrak{h}}^{s\,s}$
and $\left(-g\right)\,=\,\epsilon^{-2}$. Also, it is to be noted that
$\lim\limits_{\epsilon\,\rightarrow\,0}d^{3}y\,\tilde{\Lambda}^{t\,t}\,=\,\left(\epsilon^{6}\,d^{3}X_{L}\right)\,\left(\epsilon^{-2}\,\tilde{\Lambda}^{s\,s}\right)\,\sim\,\mathscr{O}\left(\epsilon^{4}\right)$
under the scalings (2) and (6). Therefore, under the combined temporal and
spatial coordinate scalings, the body-zone integral becomes

$\displaystyle\tilde{\mathfrak{h}}^{s\,s}\left(s,\,x^{i}\right)$
$\displaystyle=$ $\displaystyle
4\,\epsilon^{4}\,\sum_{L=1,\,2}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{\tilde{\Lambda}^{s\,s}}{\left|Z^{i}_{L}\left(s\right)\,-\,\epsilon^{2}\,X^{i}_{L}\right|}$
(68)

Now utilizing the definition of $\tilde{\Lambda}^{\mu\nu}$ from Eq. (54), at
the leading order, that is just considering the part of
$\tilde{\Lambda}^{\mu\nu}$ proportional to the classical matter energy-
momentum-stress density tensor, one obtains

$\displaystyle\tilde{\mathfrak{h}}^{s\,s}$ $\displaystyle=$ $\displaystyle
4\,\epsilon^{4}\,\sum_{L=1,\,2}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{T^{s\,s}\,\left(1\,-\,\epsilon\,f^{\prime\prime}\,\mathfrak{R}\right)}{\left|Z^{i}_{L}\,-\,\epsilon^{2}\,X^{i}_{L}\right|}$
(69)

$\epsilon$ series of the first part of the integral out of the two terms above
lead to Eq. (22), same as in GR. However, a more curious thing occurs when one
considers the second term with $\mathfrak{R}\,=\,\mathfrak{R}_{hom}$. That is
the following

$\displaystyle\tilde{\mathfrak{h}}^{s\,s}_{extra}$ $\displaystyle=$
$\displaystyle-4\,\epsilon^{5}\,f^{\prime\prime}\,\sum_{L=1,\,2}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{T^{s\,s}\,\mathfrak{R}_{hom}}{\left|Z^{i}_{L}\,-\,\epsilon^{2}\,X^{i}_{L}\right|}$
(70)

Considering now the particular form of the homogeneous Ricci density (66), one
obtains the following integral using the scaled coordinate (6)

$\displaystyle\tilde{\mathfrak{h}}^{s\,s}_{extra}$ $\displaystyle=$
$\displaystyle-4\,\epsilon^{3}\,f^{\prime\prime}\,\sum_{L=1,\,2}\,C_{L}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{T^{s\,s}}{\left|Z^{i}_{L}\left(s\right)\,-\,\epsilon^{2}\,X^{i}_{L}\right|\,\left|X^{i}_{L}\right|}\,e^{-\gamma\,\epsilon^{\frac{3}{2}}\,\left|X^{i}_{L}\right|}$
(71)

which seems to preceed the order at which the gothic metric deviation first
appears, that is $\mathscr{O}\left(\epsilon^{4}\right)$. However, under an
$\epsilon$ series about $0$, the leading order form of the above becomes the
following

$\displaystyle\tilde{\mathfrak{h}}^{s\,s}_{extra}$ $\displaystyle=$
$\displaystyle-4\,\epsilon^{3}\,f^{\prime\prime}\,\sum_{L=1,\,2}\,\frac{C_{L}}{\left|Z^{i}_{L}\right|}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{T^{s\,s}}{\left|X^{i}_{L}\right|}$
(72) $\displaystyle\tilde{\mathfrak{h}}^{s\,s}_{extra}$ $\displaystyle=$
$\displaystyle-4\,\epsilon^{3}\,f^{\prime\prime}\sum_{L=1,\,2}\,\frac{C_{L}\,N_{L}}{\left|X^{i}_{L}\right|}$
(73)

which is a Newtonian potential like term with negative one multipole moment of
the source, denoted by $N_{L}$. Therefore, the modified leading order metric
deviation $\tilde{\mathfrak{h}}^{s\,s}$ can be written in the following form

$\displaystyle\tilde{\mathfrak{h}}^{s\,s}$ $\displaystyle=$ $\displaystyle
4\,\epsilon^{4}\,\sum_{L=1,\,2}\,\left[\frac{\tilde{P}^{s}_{L}}{\left|Z^{i}_{L}\right|}\,+\,\mathscr{O}\left(\epsilon^{2}\right)\right]$
(74) $\displaystyle\tilde{P}^{s}_{L}$ $\displaystyle=$ $\displaystyle
P^{s}_{L}\,-\,\lim\limits_{\epsilon\,\rightarrow\,0}\,\frac{f^{\prime\prime}\,C_{L}}{\epsilon}\,\int_{\mathscr{B}_{L}}\,\,d^{3}X_{L}\,\frac{T^{s\,s}}{\left|X^{i}_{L}\right|}$
(75)

where $P^{s}_{L}$ was defined in Eq. (25), and is the classical definition of
quasi-local mass at the leading / Newtonian order. From astrophysical
observations, for example from solar system tests of gravity, a two-body
system obeys Newton’s law of gravity quite well, especially when their
relative velocity divided by the speed of light tends to zero. In Newton’s
law, the only parameter that arises is the classical quasi-local mass
($\frac{c^{2}\,P^{s}_{L}}{G}$ in non-geometrized units). More importantly, it
is the very first / leading order parameter that arises in a PN expansion of
GR, which has a very well defined $\epsilon\to 0$ limit, and which has been
shown to hold rather well, even in situations like a binary neutron star
system. In the current formalism, the classical qausi-local mass definition
first appears at $\mathscr{O}\left(\epsilon^{4}\right)$ for a PN expanded
$\mathfrak{h}^{s\,s}$ in GR. Noticing in Eq. (75), that $\tilde{P}^{s}_{L}$
can be interpreted as a redefined mass for the leading order
$\tilde{\mathfrak{h}}^{s\,s}$, it is important to note that while the integral
in the second term of Eq. (75) has a well defined limit for
$\epsilon\,\rightarrow\,0$, the factor multiplying it does not. In fact, due
to the $\epsilon^{-1}$ nature of the factor, the second term as a whole blows
up in the $\epsilon\,\rightarrow\,0$ limit. In order to avoid such a diverging
behavior for the classical quasi-local mass definition for the small
$\epsilon$ limit, the only choice one has is the following

$\displaystyle C_{L}$ $\displaystyle=$ $\displaystyle 0,$ (76)

which makes the homogeneous Ricci density vanish at all orders of $\epsilon$,
leading to

$\displaystyle\tilde{P}^{s}_{L}$ $\displaystyle=$ $\displaystyle P^{s}_{L}$
(77)

implying the leading order behavior of the new gothic metric deviation
$\tilde{\mathfrak{h}}^{\mu\nu}$ is the same as in GR. It is to be noted that
such a limiting procedure to Newton’s law is also used to fix the constant of
integration for the homogeneous spherically symmetric solution of GR (the
Schwarzschild mass). The vanishing of the homogeneous solution leaves one with
just the inhomogeneous/particular solution, and Eq (59) can now be solved in
an order reduced method for the Ricci density $\mathfrak{R}$ sourced by
$\tilde{\Lambda}$, as a Klein-Gordon equation with a scalar effective source
respectively. The ”non-linear terms” and the derivation of the above from Eq.
(41) has been given in Appendix A.

Substantial amount of literature is dedicated to finding non-trivial (non
Ricci flat) black hole solutions, that simultaneously solve the system of
homogeneous partial differential equations (59) and (53). See for example
constant curvature black hole space-times obtained in [47], solutions obtained
from perturbing the Schwarzschild space-time in [48], and other non-trivial
solutions in [49]. Also see [50] for spherically symmetric electro-vacuum
solutions.

Given the solution is being presented only around $R\,=\,0$, imply that at the
first approximation, only the classical matter energy-momentum tensor trace
$T_{m}$ source the Ricci scalar density $\mathfrak{R}$. The structure of the
differential equation (59), and the effective source in Eq. (60), consists of
a constant $f^{\prime\prime}\left(0\right)$ with dimension of
$\left(length\right)^{2}$, and an inverse length scale $\gamma$. Expansion
about $R\,=\,0$, and the assumption that the constant
$f^{\prime\prime}\left(0\right)$ has to be small compared to typical
$\left(\text{length}\right)^{2}$ scales of the problem (from observations
telling one that nature follows GR quite well) imply that in the RHS of Eq.
(59) the trace of the classical matter energy-momentum-stress density tensor
affect the dynamics of $\mathfrak{R}$ more strongly (that is, leading order in
the PN expansion) than terms in the curly brackets of $\tilde{\Lambda}$, in
Eq. (60), weighed by $f^{\prime\prime}\left(0\right)$. The inhomogenoeus
problem can hence be simplified into the following at the leading order

$\displaystyle\Box\,\mathfrak{R}\,+\,\frac{\gamma^{2}}{\epsilon}\,\mathfrak{R}$
$\displaystyle=$ $\displaystyle-8\,\pi\,\left(-g\right)\,\gamma^{2}\,T_{m}$
(78)

Eqs. (53) and (59) form the complete set of equations required to obtain a PN
sequence of solutions for $\mathfrak{\tilde{h}}^{\mu\nu}$ and $\mathfrak{R}$.

###  III.4 The Ricci scalar and its effect on the equations of motion

The net Ricci solution, including the particular solution, can be written as a
series about $\epsilon\,\rightarrow\,0$, and is given by (derivation in
Appendix B)

$\displaystyle\mathfrak{R}\left(t,\,x^{i}\right)$ $\displaystyle\approx$
$\displaystyle-8\,\pi\,\gamma^{2}\,\sum_{L,\,\ell,\,m}^{\ell\,=\,\infty}\,\epsilon^{2\,\ell\,+\,4}\,\frac{\left(\frac{\gamma^{2}}{\epsilon}\,-\,m^{2}\,\Omega^{2}\right)^{\frac{2\,\ell\,+\,1}{4}}}{2^{\ell\,+\,\frac{1}{2}}\,\Gamma\left(\ell\,+\,\frac{3}{2}\right)}\,\,\frac{K_{\ell\,+\,\frac{1}{2}}\left(\sqrt{\frac{\gamma^{2}}{\epsilon}\,-\,m^{2}\,\Omega^{2}}\,\left|Z^{i}_{L}\right|\right)}{\sqrt{\left|Z^{i}_{L}\right|}}$
(79)
$\displaystyle\times\,Y_{\ell\,m}\left(Z^{\theta}_{L},\,Z^{\phi}_{L}\right)\,\mathfrak{M}_{L\,\ell\,m}\,\,e^{i\,m\,\left(\phi_{0}\,-\,\Omega\,t\right)}$

where $\Gamma\left(n\right)$ is the Euler-Gamma function,
$K_{\ell\,+\,\frac{1}{2}}$ is the modified Bessel function of the second kind,
and $Y_{\ell\,m}$ are spherical harmonic functions. $\Omega$ corresponds to
the quasi-adiabatic orbital angular frequency and $\phi_{0}$ corresponds to
the phase lag between the orbital phase $\Omega\,t$ and the corresponding
phase of the Ricci scalar deviation. The phase lag is a result of the delay in
response between the movement of the sources and the corresponding Ricci
deviation generated by the accelerated motion of two ultra-compact object
about each other. The scalar multipole moments $\mathfrak{M}_{L\,\ell\,m}$ are
defined as follows

$\displaystyle\mathfrak{M}_{L\,\ell\,m}$ $\displaystyle=$ $\displaystyle
2\,\pi\,\left(-1\right)^{m}\,N_{\ell\,-m}\,\int_{\mathscr{B}_{L}}\left|X^{i}_{L}\right|^{2+\ell}\,\sin\left(X^{\theta}_{L}\right)$
$\displaystyle\times\quad P_{\ell\,-m}\left(\cos
X^{\theta}_{L}\right)\,\mathfrak{T}_{m}\left(\left|X^{i}_{L}\right|,\,X^{\theta}_{L}\right)d\left|X^{i}_{L}\right|\,dX^{\theta}_{L}$
$\displaystyle N_{\ell\,m}$ $\displaystyle=$
$\displaystyle\left(-1\right)^{m}\,\sqrt{\frac{\left(2\,\ell\,+\,1\right)}{4\,\pi}\,\frac{\left(\ell\,-\,m\right)!}{\left(\ell\,+\,m\right)!}}$
(81)

where $P_{\ell\,m}\left(X^{\theta}_{L}\right)$ are the associated Legendre
polynomials. It is to be noted that $\mathfrak{M}_{L\,0\,0}$, or the monopole
term is the usual definition of mass, as in GR, as defined in Eq. (37).

It is seen that the Ricci scalar density scales as $\epsilon^{4}$ at the
dominant order $\left(\ell\,=\,0\right)$, and its dynamics is independent of
the dynamics of the modified gothic metric $\mathfrak{\tilde{h}}^{\mu\nu}$. It
is however, for the sake of obtaining observational consequences of the
current theory, more interesting to look at the opposite, that is the effect
of the Ricci scalar deviation on the dynamics of
$\mathfrak{\tilde{h}}^{\mu\nu}$, and consequently, its gauge invariant
contribution to the modified equations of motion of the binary system.

The equations of motion can be found by using the modified conservation law
(55) alongwith the definitions of (II.2) and (36) with $\Lambda^{\mu\nu}$
replaced by $\tilde{\Lambda}^{\mu\nu}$. Since the equations of motion involve
a vector surface integral over the body zone boundary, on which
$T_{m}^{\mu\nu}$ vanishes (the classical matter source being a compact one),
and consequently so does the effect of the first term in the RHS of Eq. (56)
in the equations of motion. A vector surface integral over a sphere also
vanishes if the integrand is spherically symmetric, which
${}_{(4)}\mathfrak{R}$ is, as seen in (116). Therefore, the terms that might
contribute to a deviation from GR at the leading order are the space-space
components of the effective energy-momentum tensor $t^{\mu\nu}_{eff}$, given
by

$\displaystyle
t^{i\,j}_{eff}=\frac{\left(-g\right)^{-2}}{192\,\pi\,\gamma^{4}}\,\left[2\,\left(\eta^{i\,k}\,\eta^{j\,l}\,+\,\eta^{i\,l}\,\eta^{j\,k}\,-\,\eta^{i\,j}\,\eta^{k\,l}\right)\,\mathfrak{R}\,\partial_{k}\mathfrak{R}\,\partial_{l}\tilde{\mathfrak{h}}\right.$
$\displaystyle\left.\,+\,\left(\eta^{i\,k}\,\eta^{j\,l}\,-\,\frac{1}{2}\,\eta^{i\,j}\,\eta^{k\,l}\right)\,\left(4\,\partial_{k}\mathfrak{R}\,\partial_{l}\mathfrak{R}\,+\,\mathfrak{R}^{2}\,\partial_{k}\tilde{\mathfrak{h}}\,\partial_{l}\tilde{\mathfrak{h}}\right)\right]$
(82)

owing to the fact that at the leading order, both the metric density deviation
trace $\tilde{\mathfrak{h}}$ and the Ricci deviation density $\mathfrak{R}$
are time independent, and $\left|\eta^{s\,s}\right|\,\sim\,\epsilon^{2}$.

One notices that the leading and sub-leading order Ricci density scale as
$\mathscr{O}\left(\epsilon^{4}\right)$, and correspondingly, the deviation
from GR in the equations of motion of the binary appears at
$\mathscr{O}\left(\epsilon^{12}\right)$, which leads to the modification of
the 6 PN equations of motion, which are then given by the following

$\displaystyle M_{1}\,\frac{dv_{1}^{i}}{ds}$ $\displaystyle=$ $\displaystyle
F^{i}_{Newtonian}\,+\,\epsilon^{2}\,F^{i}_{1PN,GR}\,+\,\epsilon^{4}\,F^{i}_{2PN,GR}\,+\,\epsilon^{5}\,F^{i}_{2.5PN,GR}$
(83)
$\displaystyle\,+\,F^{i}_{3PN,GR}\,+\,...\,+\,\epsilon^{12}\,\left(F^{i}_{6PN,GR}\,-\,\frac{32\,\pi\,M_{1}\,M_{2}^{3}}{9\,r_{1\,2}^{4}}\,\,e^{-\frac{2\,\gamma\,r_{1\,2}}{\sqrt{\epsilon}}}\,n^{i}\right)$
$\displaystyle\,+\,\mathscr{O}\left(\epsilon^{13}\right)$

The derivation of the above is given in Appendix B. $F_{N}^{i}$ and
$F_{kPN,GR}^{i}$ are the attractive Newtonian and
$k^{th}\,\left(k\,=\,1,...,6\right)$ PN forces in GR, respectively, on the
first body zone.

The leading modification at 6 PN to the equations of motion is seen to be
repulsive at very short scales, and for very high velocities and compactness,
given the convention chosen in this paper regarding the direction of the unit
normal $n^{i}$, which points from the first body zone towards the COM of the
second source. The leading modification is dependent on a length scale that
acts like a new universal constant (if nature follows such a quadratic theory)
and is $\sim\,\gamma^{-1}$. It is interesting to note that in the 6 PN extra
force term, one can cancel $M_{1}$ from both sides leading to only $M_{2}^{3}$
remaining as a product in the entire RHS of Eq. (83). This establishes that
the equivalence/effacing nature of two ultra-compact sources in a binary orbit
around each other till the 6 PN order in polynomial $f\left(R\right)$ theories
of gravity with $f^{\prime\prime}\left(0\right)\,<\,0$ remains almost the same
as in GR.

##  IV Discussions and Conclusions

In this article the PN formalism of GR was reviewed using the surface integral
approach [24], using the specially scaled co-moving coordinate system for the
internal problem [14]. Using the strong field point particle mechanism
developed for the PN expansion of GR in [14], a quadratic $f\left(R\right)$
model of gravity was PN expanded in the conformal-harmonic gauge. In doing so,
it was clearly established that no deviation from GR occurs in the quadratic
$f\left(R\right)$ theory till the sixth post Newtonian order, thus
automatically also taking care of the effacing nature of sources in a binary
orbit.

A scaled co-moving coordinate system, and the Einstein-Infeld approach [24]
was used in order to solve for the two body equations of motion in
$f\left(R\right)$ theories of gravity. A conformal scaling of the old
definition of the gothic metric of GR was found to be more suited to get
possible observational consequences out of a theory that is highly non-linear
and fourth derivative of the metric tensor. Under the conformal scaling, the
field equations resemble the familiar equations of the Landau-Lifshitz
formulation of GR, with an extra ’effective’ source term that is the
manifestation of the extra scalar degree of freedom in $f\left(R\right)$
theories.

More notably, a dynamical scalar variable for solving the combined
tensor+scalar system of second order partial differential equations was found:
the Ricci scalar density $\mathfrak{R}$. This particular redefinition frames
the tensor differential equations and the scalar Klein-Gordon equation in a
similar footing; that is expanded in a series of
$\frac{f^{\prime\prime}}{\sqrt{-g}}$, that is, $\epsilon\,f^{\prime\prime}$.
Coupling the scalar variable, the Ricci scalar, with the volume element
$\sqrt{-g}$ in the harmonic gauge of Eq. (51) leads to a system of
differential equations that can be solved perturbatively in terms of the PN
parameter $\epsilon$, that smoothly evolves from zero to unity. The
redefinition allows one to ignore the Klein-Gordon homogeneous solution; and
to sieve through non-linear terms of $\tilde{\mathfrak{h}}^{\mu\nu}$ and
$\mathfrak{R}$ (and their first derivatives) in the effective energy-momentum
tensor $t^{\mu\nu}_{eff}$, and obtain terms that are leading with regards to
the parameter $\epsilon$, as well as group them in terms of sub-leadingness.

It was shown that the homogeneous Ricci scalar density has to vanish if the
definition of mass at the Newtonian level is well defined. This then leads to
the inhomogeneous Ricci scalar density solution appearing at
$\mathscr{O}\left(\epsilon^{4}\right)$, and sourcing the metric density tensor
deviation $\tilde{\mathfrak{h}}^{\mu\nu}$ through the effective source term
$t^{\mu\nu}_{eff}$ of $f\left(R\right)$ theories of gravity.

In this study, the change in the equations of motion of two extended (but
compact objects) from GR in the quadratic Ricci model till the $6^{th}$ PN
order were obtained. The equations of motion till 5.5 PN order remain
unchanged compared to GR. A modification of the equations of motion from the
GR form, through the appearance of a repulsive Yukawa like force at length
scales comparable to $\gamma^{-1}$ at the $6^{th}$ PN order was observed.
However, the modified force is time independent, conservative, and dependent
only on the usual definition of mass & the orbital separation. A universal
constant $f^{\prime\prime}\left(0\right)$, other than the gravitational
constant (put to unity in Geometrized units), was defined, similar to other
works in the literature. The constant has dimensions of
$\left(length\right)^{2}$, and a corresponding inverse length scale $\gamma$.

The above analysis also brings to light a possibly new class of solutions of
$f\left(R\right)$ theories of gravity, and shows that the appearance of
structure in the equations of motion in $f\left(R\right)$ theories are exactly
the same as in GR, to a very high PN accuracy. At the leading order of
deviation from GR, no extra structure comes into the equations of motion due a
quadratic change in the action. The appearance of structure at high PN orders
till 6 PN is still solely a GR effect. This is the strong equivalence
principle of $f\left(R\right)$ theories in action. The current study
establishes the proposition that both the weak and strong equivalence
principles hold for $f\left(R\right)$ theories of gravity, and tests for the
validity of GR through tests of the equivalence principle will not at all
distinguish between the two theories. Hence, the equivalence, and
consequently, the effacing principle is not unique to GR.

## Acknowledgement

The author would like to thank Bruce Allen and Badri Krishnan for a
significant number of ideas and clarifications throughout the project. The
author thanks Pierre Mourier for important comments on earlier drafts. The
author would like to thank Yousuke Itoh, and the CBC group of MPI for
Gravitational Physics, especially Andrey Shoom, for important comments and
questions. The author would like to thank Chowhan Tanmay Tushar and Amanda
Fajardo Garcia for some key inputs to the calculations of the Landau-Lifshitz
formalism of GR. Finally, the author would also like to thank the anonymous
reviewer for suggesting important changes to the manuscript throughout the
review process.

The technical computing system Mathematica, and some of the symbolic tensor
manipulation packages supported by it, namely $xAct$ [51] was used for most of
the calculations in the study.

##  Appendix A Taming the scalar wave equation for the extra massive degree of
freedom of $f\left(R\right)$

The scalar wave dynamics of Eq. (41) is non-linear, as seen in terms like
$f^{\prime}R$. By using the definition of the Laplace-Beltrami operator, one
may write (53) as

$\displaystyle\frac{3}{\sqrt{-g}}\,\left(\sqrt{-g}\,g^{\mu\nu}\,f^{\prime}_{\,,\,\nu}\right)_{\,,\,\mu}\,+\,\left(f^{\prime}\,R\,-\,2\,f\right)$
$\displaystyle=$ $\displaystyle 8\,\pi\,T_{m}$

Noticing the old gothic metric definition $\left(\sqrt{-g}\,g^{\mu\nu}\right)$
in the first term, one can conformally transform the old definition to the new
definition using Eq. (44) and perform the perturbation in Eq. (50), after
which the above equations become

$\displaystyle-\frac{3}{\sqrt{-g}\,f^{\prime}}\,\left[\Box
f^{\prime}\,-\,\mathfrak{\tilde{h}}^{\mu\nu}\,f^{\prime}_{\,,\,\mu\nu}\,-\,\frac{1}{f^{\prime}}\,\left\\{\left(f^{\prime}\right)^{\,,\,\mu}\,f^{\prime}_{\,,\,\mu}\right.\right.$
$\displaystyle\left.\left.+\,\mathfrak{\tilde{h}}^{\mu\nu}\,f^{\prime}_{\,,\,\mu}f^{\prime}_{\,,\,\nu}\right\\}\right]\,+\,\left(f^{\prime}\,R\,-\,2\,f\right)=8\,\pi\,T_{m}$
(85)

The second and the fourth term inside the square brackets can be transformed,
by chain rule, into total derivative terms, and terms where the Lorenz gauge
conditions can be imposed. Noting that boundary dependent terms can be
discarded to obtain gauge invariant objects, one obtains

$\displaystyle\Box
f^{\prime}\,+\,\sqrt{-g}\left[\frac{1}{3}\,\left(f^{\prime}\right)^{2}\,R\,-\,\frac{2}{3}\,f\,f^{\prime}\,-\,\frac{8\,\pi}{3}\,f^{\prime}\,T_{m}\right]$
$\displaystyle=$ $\displaystyle 0$

Now the particular forms for $f\left(R\right)$ and $f^{\prime}\left(R\right)$
can be imposed, according to Eqs. (48) and (49), as well as the dynamical
variable definition for the Ricci scalar density $\mathfrak{R}$ in Eq. (52),
which leads to

$\displaystyle\Box\,\mathfrak{R}\,+\,\sqrt{-g}\,\gamma^{2}\,\mathfrak{R}=-8\,\pi\,\left(-g\right)\,\gamma^{2}\,\tilde{\Lambda}$
(87)
$\displaystyle\tilde{\Lambda}=\,T_{m}\,+\,f^{\prime\prime}\left(0\right)\,\left\\{\frac{T_{m}\,\mathfrak{R}}{3\,\sqrt{-g}}+\,\frac{1}{8\,\pi\,\left(-g\right)}\left(\frac{\mathfrak{R}^{2}}{3}\,-\,\partial_{\mu}\mathfrak{R}\,\partial^{\mu}\tilde{\mathfrak{h}}\,-\,\frac{1}{2}\,\mathfrak{R}\,\partial_{\mu}\tilde{\mathfrak{h}}^{\alpha\beta}\,\partial^{\mu}\tilde{\mathfrak{h}}_{\alpha\beta}\right.\right.$
$\displaystyle\left.\left.\,-\,\frac{1}{4}\,\mathfrak{R}\,\partial_{\mu}\tilde{\mathfrak{h}}\,\partial^{\mu}\tilde{\mathfrak{h}}\,-\,\frac{1}{2}\,\mathfrak{R}\,\Box\tilde{\mathfrak{h}}\right)\right\\}\,+\,\mathscr{O}\left(\left[f^{\prime\prime}\left(0\right)\right]^{2}\right)$
(88)
$\displaystyle\gamma^{2}\equiv-\frac{1}{3\,f^{\prime\prime}\left(0\right)}$
(89) $\displaystyle\,T_{m}\,=\,\eta_{\mu\nu}\,T_{m}^{\mu\nu}$ (90)

where the argument of $\mathscr{O}$ in the above have been made dimensionless.
Since in the scaled time, inverse of the square root of the determinant of the
metric scale as $\epsilon$, corrections to the effective energy-momentum
scalar $\tilde{\Lambda}$ from $T_{m}$ at the next-to-leading order will be of
$\mathscr{O}\left[\epsilon^{2}\,\frac{f^{\prime\prime}\left(0\right)}{\left(\text{length}\right)^{2}}\right]$,
with $\left(\text{length}\right)$ being typical length scales of the binary
problem.

##  Appendix B Expanding the solution of the order reduced Klein-Gordon
equation in orders of $\epsilon$

The net solution to the Klein-Gordon problem (78) with the boundary condition
choice of no incoming scalar radiation from past null Minkowskian infinity is
then given as follows

$\displaystyle\mathfrak{R}\left(x^{\mu}\right)$ $\displaystyle=$
$\displaystyle-8\,\pi\,\gamma^{2}\,\epsilon^{-2}\,\int\,d^{4}y\,G_{\tilde{\gamma}}\left(x^{\mu},\,y^{\mu}\right)\,T_{m}\left(y^{\mu}\right)$
(91) $\displaystyle\tilde{\gamma}^{2}$ $\displaystyle\equiv$
$\displaystyle\frac{\gamma^{2}}{\epsilon}$ (92)

where the factor of $\epsilon^{-2}$ arises because the metric determinant
$\left(-g\right)$ multiplying the RHS of Eq. (78), under the transformation
(2), scale as $\epsilon^{-2}$, as seen after Eq. (58).
$G_{\tilde{\gamma}}\left(x^{\mu},\,y^{\mu}\right)$ is the retarded Green’s
function of the Klein-Gordon equation, as was given in [40] with the
$\left(+\,,\,-\,,\,-\,,\,-\right)$ metric signature as

$\displaystyle G_{\tilde{\gamma}}\left(t,\,q;\,x^{i},\,y^{i}\right)$
$\displaystyle=$
$\displaystyle\int_{-\infty,\,\tilde{\gamma}}^{-\tilde{\gamma},\,\infty}\frac{d\,\omega}{2\,\pi}\,e^{-i\,\omega\left(t-q\right)}\,\frac{e^{i\,\sqrt{\omega^{2}-\tilde{\gamma}^{2}}\,\left|x^{i}-y^{i}\right|}}{4\,\pi\,\left|x^{i}-y^{i}\right|}$
$\displaystyle+\,\int_{-\tilde{\gamma}}^{\tilde{\gamma}}\,\frac{d\omega}{2\,\pi}\,e^{-i\,\omega\,\left(t-q\right)}\,\frac{e^{-\,\sqrt{\tilde{\gamma}^{2}-\omega^{2}}\,\left|x^{i}-y^{i}\right|}}{4\,\pi\,\left|x^{i}-y^{i}\right|}$

The notation $\int_{-\infty,\,\tilde{\gamma}}^{-\tilde{\gamma},\,\infty}$
involve two integrals, one from $-\infty$ to $-\tilde{\gamma}$, and the other
from $\tilde{\gamma}$ to $\infty$. It is to be noted that the three integrals
whose domains encompass all of $\omega$ space must be evaluated simultaneously
in order for the solution to converge. $T_{m}\left(q,\,y^{i}\right)$, the
trace of the classical matter energy-momentum tensor, which being a scalar,
does not transform under any of the coordinate scalings.

It is to be noted that the PN parameter $\epsilon$ appears explicitly in the
Green’s function, both as an argument of the exponential function, and as one
(or both) of the limits of the two $\omega$ space integrals at hand.
Considering the first integral, that is

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\int_{-\infty,\,\tilde{\gamma}}^{-\tilde{\gamma},\,\infty}\frac{d\,\omega}{2\,\pi}\,e^{-i\,\omega\left(t-q\right)}\,\frac{e^{i\,\sqrt{\omega^{2}-\tilde{\gamma}^{2}}\,\left|x^{i}-y^{i}\right|}}{4\,\pi\,\left|x^{i}-y^{i}\right|}$
(94)

In the Newtonian limit, or $\epsilon\to 0$, $\tilde{\gamma}$ diverges. But
since it appears in the argument of the oscillatory function, that is
$e^{i\,\sqrt{\omega^{2}-\tilde{\gamma}^{2}}\,\left|x^{i}-y^{i}\right|}$, the
Green’s function itself does not diverge and remains finite in
$\left[-1,1\right]$ as long as $\omega^{2}>\tilde{\gamma}^{2}$. Although it
becomes a highly oscillating function. At the same time the limits of the
integral are replaced with $\pm\tilde{\gamma}\to\pm\infty$, effectively
shrinking the domain of integration to a vanishing value (implying that the
lower and upper limits of the two integrals coincide for $\epsilon\to 0$).
This leads to a vanishing value of $I_{1}$ for $\epsilon\to 0$. The second
integral, that is when $\tilde{\gamma}^{2}>\omega^{2}$, which is the following

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\int_{-\tilde{\gamma}}^{\tilde{\gamma}}\,\frac{d\omega}{2\,\pi}\,e^{-i\,\omega\,\left(t-q\right)}\,\frac{e^{-\,\sqrt{\tilde{\gamma}^{2}-\omega^{2}}\,\left|x^{i}-y^{i}\right|}}{4\,\pi\,\left|x^{i}-y^{i}\right|},$
(95)

The limit $\epsilon\to 0$ leads to the domain of integration covering the
whole of Real number space. However, the argument in the decaying function
$e^{-\,\sqrt{\tilde{\gamma}^{2}-\omega^{2}}\,\left|x^{i}-y^{i}\right|}$ blows
up leading to the exponential vanishing, and hence leading to $I_{2}\to 0$ for
$\epsilon\to 0$.

Due to the the coordinate scalings in Eq. (6), the infinitesimal 4-volume
element $d^{4}y$ transform from the asymptotic observer’s frame to either of
the body zones $\mathscr{B}_{L}$ in the following manner

$\displaystyle d^{4}y$ $\displaystyle\equiv$ $\displaystyle
dt\,\wedge\,d^{3}y$ (96) $\displaystyle\rightarrow$
$\displaystyle\epsilon^{-1}\,ds\,\wedge\,\epsilon^{6}\,\,d^{3}X_{L}$ (97)

where $\wedge$ denotes the wedge product between 1-form dt and 3-form
$d^{3}y$. In the current operational context $\wedge$ is effectively scalar
multiplication.

Substituting Eq. (B) in Eq. (91), and transforming into the body zone
coordinates by substituting Eq. (97) as the infinitesimal covariant volume
element in Eq. (91), one obtains the following integral for the particular
solution $\mathfrak{R}_{part}$ for the inhomogeneous Klein-Gordon equation

$\displaystyle\mathfrak{R}_{part}\left(t,\,x^{i}\right)$ $\displaystyle=$
$\displaystyle-\epsilon^{3}\,\,8\,\pi\,\gamma^{2}\,\int_{-\infty}^{\infty}\,ds^{\prime}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\,\int_{-\infty,\,\tilde{\gamma}}^{-\tilde{\gamma},\,\infty}\,\frac{d\omega}{2\,\pi}\,e^{-i\,\omega\,\left(t-s^{\prime}/\epsilon\right)}$
$\displaystyle\times\quad\,\frac{e^{i\,\sqrt{\omega^{2}-\tilde{\gamma}^{2}}\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}}{4\,\pi\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}\,\left[\epsilon^{4}\,T_{m}\left(s^{\prime}/\epsilon,\,X_{L}^{i}\right)\right]$
$\displaystyle-\,\epsilon^{3}\,\quad\,8\,\pi\,\gamma^{2}\,\int_{-\infty}^{\infty}\,ds^{\prime}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\int_{-\tilde{\gamma}}^{\tilde{\gamma}}\,\frac{d\omega}{2\,\pi}\,e^{-i\,\omega\,\left(t-s^{\prime}/\epsilon\right)}$
$\displaystyle\times\quad\,\frac{e^{-\sqrt{\tilde{\gamma}^{2}-\omega^{2}}\,\left|Z_{L}^{i}-\epsilon^{2}\,X^{i}_{L}\right|}}{4\,\pi\,\left|Z_{L}^{i}-\epsilon^{2}\,X^{i}_{L}\right|}\,\left[\epsilon^{4}\,T_{m}\left(s^{\prime}/\epsilon,\,X^{i}_{L}\right)\right]$

In the body zones that surround each of the objects, the energy-momentum
tensor of classical matter
$\left(T_{\mathscr{B}_{L}}\,=\,\lim\limits_{\epsilon\,\rightarrow\,0}\epsilon^{4}\,T_{m}\right)$
will be assumed to have a quasi-stationary/adiabatic initial condition for
solving the relaxed system of equations, in the absence of any other time
dependent driving force inside the compact object, and is given by

$\displaystyle T_{\mathscr{B}_{L}}\left(X^{i}\right)$ $\displaystyle=$
$\displaystyle\sum_{n=-\infty}^{\infty}\,\mathfrak{T}_{n}\left(\left|X^{i}_{L}\right|,X^{\theta}_{L}\right)\,e^{i\,n\,\mathscr{X}}$
(99) $\displaystyle\mathscr{X}\left(s/\epsilon,\,X_{L}^{\phi}\right)$
$\displaystyle\equiv$
$\displaystyle\mathscr{X}\,=\,X^{\phi}_{L}\,-\,\Omega\,s/\epsilon\,+\,\phi_{0}$
(100)

where an axial+time symmetry was assumed for the initial condition, such that
the time dependence of the source body in the body zone coordinate system
repeats after every $T\,=\,\frac{2\,\pi}{\Omega}$, with a constant phase
parameter $\phi_{0}$. The weighing factors $\mathfrak{T}_{n}$ are coefficients
in the series expansion of the trace $T_{m}\left(\text{or
}T_{\mathscr{B_{L}}}\right)$ using stationary functions
$e^{i\,n\,\mathscr{X}}$, and are functions of the radial and azimuthal
coordinates. Eq. (99) physically implies that the variations in the energy-
momentum tensor sourcing $\mathfrak{R}$ in the body zone of the first object
is purely generated by the effect of the motion of the second object around
it. The body zone coordinates $\left(X^{i}_{L}\right)$ were defined in Eq.
(6), in which $\left|X^{i}_{L}\right|$ is the distance from the center of mass
of the body $L$ to any point in the body zone coordinates, as viewed in the
respective body zones. The choice in Eq. (99) has a simplifying effect on the
subsequent calculations and is justified by the adiabatic and stationarity in
the co-moving frame approximations, as found in the literature on PN
expansions.

Assuming quasi-periodicity for the source tensor, as in Eq. (99), one can
substitute it into Eq. (B), to obtain

$\displaystyle\mathfrak{R}_{part}\left(t,\,x^{i}\right)$ $\displaystyle=$
$\displaystyle-\epsilon^{3}\,\,\frac{8\,\pi\,\gamma^{2}}{2\,\pi}\,\sum_{n=-\infty}^{\infty}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\left(\int_{-\infty,\,\tilde{\gamma}}^{-\tilde{\gamma},\,\infty}\,d\omega\,e^{-i\,\omega\,t}\right)\,e^{i\,n\,\left(X^{\phi}_{L}+\phi_{0}\right)}$
(101)
$\displaystyle\times\quad\,\frac{e^{i\,\sqrt{\omega^{2}-\tilde{\gamma}^{2}}\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}}{4\,\pi\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}\,\mathfrak{T}_{n}\left(\left|X^{i}_{L}\right|,X^{\theta}_{L}\right)\,\int_{-\infty}^{\infty}\,ds^{\prime}\,e^{i\,\left(\omega-n\,\Omega\right)\,s^{\prime}/\epsilon}$
$\displaystyle\,-\,\epsilon^{3}\quad\,\frac{8\,\pi\,\gamma^{2}}{2\,\pi}\,\sum_{n=-\infty}^{\infty}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\left(\int_{-\tilde{\gamma}}^{\tilde{\gamma}}\,d\omega\,e^{-i\,\omega\,t}\right)\,e^{i\,n\,\left(X^{\phi}_{L}+\phi_{0}\right)}$
$\displaystyle\times\quad\,\frac{e^{-\,\sqrt{\tilde{\gamma}^{2}-\omega^{2}}\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}}{4\,\pi\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}\,\mathfrak{T}_{n}\left(\left|X^{i}_{L}\right|,X^{\theta}_{L}\right)\,\int_{-\infty}^{\infty}\,ds^{\prime}\,e^{i\,\left(\omega-n\,\Omega\right)\,s^{\prime}/\epsilon}$

Transforming

$\displaystyle s^{\prime}\rightarrow\epsilon\,s$ (102)

and using the integral representation of the delta function,

$\displaystyle\int_{-\infty}^{\infty}\,e^{i\,k\,\left(x-x^{\prime}\right)}\,dk$
$\displaystyle=$ $\displaystyle 2\,\pi\,\delta\left(x-x^{\prime}\right),$
(103)

one obtains

$\displaystyle\mathfrak{R}_{part}\left(t,\,x^{i}\right)$ $\displaystyle=$
$\displaystyle-\epsilon^{4}\,\,8\,\pi\,\gamma^{2}\,\sum_{n=-\infty}^{\infty}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{e^{i\,\sqrt{n^{2}\,\Omega^{2}-\tilde{\gamma}^{2}}\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}}{4\,\pi\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}$
(105)
$\displaystyle\times\quad\,\mathfrak{T}_{n}\left(\left|X^{i}_{L}\right|,X^{\theta}_{L}\right)\,\times\,e^{-i\,n\,\left(\Omega\,t\,-\,X^{\phi}_{L}\,-\,\phi_{0}\right)}$
$\displaystyle-\epsilon^{4}\,\,8\,\pi\,\gamma^{2}\,\sum_{n=-\infty}^{\infty}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{e^{-\,\sqrt{\tilde{\gamma}^{2}-n^{2}\,\Omega^{2}}\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}}{4\,\pi\,\left|Z^{i}_{L}-\epsilon^{2}\,X^{i}_{L}\right|}$
$\displaystyle\times\quad\,\mathfrak{T}_{n}\left(\left|X^{i}_{L}\right|,X^{\theta}_{L}\right)\,e^{-i\,n\,\left(\Omega\,t\,-\,X^{\phi}_{L}\,-\,\phi_{0}\right)}$

The Green’s function (B) can be written as an infinite sum of spherical
harmonic functions, that are weighed by functions of the radial coordinates
$\left|Z^{i}_{L}\right|$ and $\left|X^{i}_{L}\right|$, in the following way
for $\left|Z^{i}_{L}\right|\,>\,\left|X^{i}_{L}\right|$

$\displaystyle\frac{e^{i\,k\,\left|X^{i}\,-\,Z^{i}\right|}}{4\,\pi\,\left|X^{i}\,-\,Z^{i}\right|}$
$\displaystyle=$ $\displaystyle
i\,k\,\sum_{\ell,\,m}j_{\ell}\left(k\,\left|X^{i}\right|\right)\,h_{\ell}^{(1)}\left(k\,\left|Z^{i}\right|\right)\,Y^{*}_{\ell\,m}\left(X^{\theta},\,X^{\phi}\right)\,Y_{\ell\,m}\left(Z^{\theta},\,Z^{\phi}\right)$
$\displaystyle\frac{e^{-\,k\,\left|X^{i}\,-\,Z^{i}\right|}}{4\,\pi\,\left|X^{i}\,-\,Z^{i}\right|}$
$\displaystyle=$
$\displaystyle\sum_{\ell,\,m}\,\frac{I_{\ell+\frac{1}{2}}\left(k\,\left|X^{i}\right|\right)\,K_{\ell+\frac{1}{2}}\,\left(k\,\left|Z^{i}\right|\right)}{\sqrt{\left|X^{i}\right|\,\left|Z^{i}\right|}}\,Y^{*}_{\ell\,m}\left(X^{\theta},\,X^{\phi}\right)\,Y_{\ell\,m}\left(Z^{\theta},\,Z^{\phi}\right)$
$\displaystyle Z^{i}_{L}$ $\displaystyle\equiv$
$\displaystyle\frac{x^{i}\,-\,z^{i}_{L}\left(s\right)}{\epsilon^{2}}\,;\quad\,L\,=\,1,\,2.$
(107)

The various functions appearing above are as follows

  * •

$j_{\ell}$: Spherical Bessel function of first kind.

  * •

$h_{\ell}^{(1)}$: Spherical Bessel function of third kind.

  * •

$I_{\ell\,+\,\frac{1}{2}}$: Modified Bessel function of the first kind.

  * •

$K_{\ell\,+\,\frac{1}{2}}$: Modified Bessel function of the second kind.

  * •

$Y_{\ell\,m},\,Y^{*}_{\ell\,m}$: Spherical harmonic functions and their
complex conjugates.

After using Eqs. (LABEL:greenexpupper) and (107) in the above, one obtains for
$\left|Z^{i}_{L}\right|\,>\,\left|X^{i}_{L}\right|$

$\displaystyle\mathfrak{R}_{part}\left(t,\,x^{i}\right)$ $\displaystyle=$
$\displaystyle-\epsilon^{4}\,\,8\,\pi\,i\,\gamma^{2}\,\sum_{n,\,\ell,\,m}^{\left|n\right|\,>\,\left\lfloor\frac{\tilde{\gamma}}{\Omega}\right\rfloor}\,\sqrt{n^{2}\,\Omega^{2}-\tilde{\gamma}^{2}}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,j_{\ell}\left(\epsilon^{2}\,\sqrt{n^{2}\,\Omega^{2}-\tilde{\gamma}^{2}}\,\left|X_{L}^{i}\right|\right)$
$\displaystyle\times\,\,h_{\ell}^{(1)}\left(\sqrt{n^{2}\,\Omega^{2}-\tilde{\gamma}^{2}}\,\left|Z^{i}\right|\right)\quad\,Y^{*}_{\ell\,m}\left(X_{L}^{\theta},\,X_{L}^{\phi}\right)\,Y_{\ell\,m}\left(Z^{\theta},\,Z^{\phi}\right)\,\mathfrak{T}_{n}\left(\left|X_{L}^{i}\right|,\,X_{L}^{\theta}\right)\,e^{i\,n\,\left(X_{L}^{\phi}+\phi_{0}-\Omega\,t\right)}$
$\displaystyle-\,\epsilon^{3}\,\quad\,\,8\,\pi\,\gamma^{2}\,\sum_{n,\,\ell,\,m}^{\left|n\right|\,<\,\left\lfloor\frac{\tilde{\gamma}}{\Omega}\right\rfloor}\,\int_{\mathscr{B}_{L}}\,d^{3}X_{L}\,\frac{I_{\ell+\frac{1}{2}}\left(\epsilon^{2}\,\sqrt{\tilde{\gamma}^{2}-n^{2}\,\Omega^{2}}\,\left|X_{L}^{i}\right|\right)}{\sqrt{\left|X_{L}^{i}\right|}}\,\frac{K_{\ell+\frac{1}{2}}\left(\sqrt{\tilde{\gamma}^{2}-n^{2}\,\Omega^{2}}\,\left|Z^{i}\right|\right)}{\sqrt{\left|Z^{i}\right|}}$
$\displaystyle\times\quad\,Y^{*}_{\ell\,m}\left(X_{L}^{\theta},\,X_{L}^{\phi}\right)\,Y_{\ell\,m}\left(Z^{\theta},\,Z^{\phi}\right)\,\mathfrak{T}_{n}\left(\left|X_{L}^{i}\right|,\,X_{L}^{\theta}\right)\,e^{i\,n\,\left(X_{L}^{\phi}+\phi_{0}-\Omega\,t\right)}$

one can substitute the following properties of spherical harmonic functions in
Eq. (B)

$\displaystyle Y^{*}_{\ell\,m}=\left(-1\right)^{m}\,Y_{\ell\,-m}$ (109)
$\displaystyle
Y_{\ell\,-m}=N_{\ell\,-m}\,P_{\ell\,-m}\,\left(\cos\,X^{\theta}_{L}\right)\,e^{-i\,m\,X^{\phi}}$
(110)
$\displaystyle\int_{0}^{2\,\pi}\,dX^{\phi}_{L}\,Y^{*}_{\ell\,m}\,e^{i\,n\,X^{\phi}}=2\,\pi\,\left(-1\right)^{m}\,N_{\ell\,-m}\,P_{\ell\,-m}\left(X^{\theta}_{L}\right)\,\delta_{m\,n}$
(111) $\displaystyle
N_{\ell\,m}=\left(-1\right)^{m}\,\sqrt{\frac{\left(2\,\ell\,+\,1\right)}{4\,\pi}\,\frac{\left(\ell\,-\,m\right)!}{\left(\ell\,+\,m\right)!}}$
(112)

where $P_{\ell\,m}\left(X^{\theta}_{L}\right)$ are the associated Legendre
polynomials and $\delta_{m\,n}$ is the Kronecker delta distribution. Since $m$
takes values between $-\ell$ to $\ell$, application of the Kronecker delta
leads to the summation on $\ell$ in the first term of Eq. (B) going from the
lower integral part of $\frac{\tilde{\gamma}}{\Omega}$ or
$\left\lfloor\frac{\tilde{\gamma}}{\Omega}\right\rfloor$ to $\infty$, whereas
the second terms’ summation range takes $\ell$ from zero to
$\left\lfloor\frac{\tilde{\gamma}}{\Omega}\right\rfloor$. The contributions to
the particular solution from above the cutoff energy scale $\tilde{\gamma}$ is
highly suppressed (especially at lower PN orders) because of the properties of
$j_{\ell}$ and $h^{(1)}_{\ell}$ for asymptotically large values of $\ell$, and
terms in $\mathfrak{R}$ series containing the same can be approximated to zero
at the current orders of calculation. Hence one can safely set the first
summation in RHS of Eq. (B) to zero and consider only the second summation.
Hence it is to be noted that the factors that depend on $\epsilon$ in the
above are only in the arguments of the Bessel function
$I_{\ell\,+\,\frac{1}{2}}$.

Hence the only relevant factors that explicitly contain $\epsilon$ in both the
sums in Eq. (B), can be expanded about $\epsilon\,\rightarrow\,0$ as follows

$\displaystyle
I_{\ell+\frac{1}{2}}\,\left(\epsilon^{2}\,\Lambda\,\left|X_{L}\right|\right)$
$\displaystyle\approx$
$\displaystyle\epsilon^{2\,\ell\,+\,1}\,\frac{\Lambda^{\frac{2\,\ell\,+\,1}{2}}}{2^{\ell\,+\,\frac{1}{2}}\,\Gamma\left(\ell\,+\,\frac{3}{2}\right)}\,\left|X_{L}\right|^{\ell\,+\,\frac{1}{2}}$

The approximation of (B), when substituted in for the relevant factors of the
RHS of Eq. (B), leads to the following polynomial series of $\epsilon$ about
$\epsilon\,\rightarrow\,0$

$\displaystyle\mathfrak{R}_{part}\left(t,\,x^{i}\right)\approx-8\,\pi\,\gamma^{2}\,\sum_{L,\,\ell,\,m}^{\ell\,=\,\infty}\,\epsilon^{2\,\ell\,+\,4}\,\frac{\left(\tilde{\gamma}^{2}\,-\,m^{2}\,\Omega^{2}\right)^{\frac{2\,\ell\,+\,1}{4}}}{2^{\ell\,+\,\frac{1}{2}}\,\tilde{\gamma}\left(\ell\,+\,\frac{3}{2}\right)}$
$\displaystyle\times\,e^{i\,m\,\left(\phi_{0}\,-\,\Omega\,t\right)}\,\,\frac{K_{\ell\,+\,\frac{1}{2}}\left(\sqrt{\tilde{\gamma}^{2}\,-\,m^{2}\,\Omega^{2}}\,\left|Z^{i}_{L}\right|\right)}{\sqrt{\left|Z^{i}_{L}\right|}}\,Y_{\ell\,m}\left(Z^{\theta}_{L},\,Z^{\phi}_{L}\right)$
$\displaystyle\times\,\mathfrak{M}_{L\,\ell\,m}$ (114)
$\displaystyle\mathfrak{M}_{L\,\ell\,m}=2\,\pi\,\left(-1\right)^{m}\,N_{\ell\,-m}\,\int_{\mathscr{B}_{L}}\left|X^{i}_{L}\right|^{2+\ell}\,\sin\left(X^{\theta}_{L}\right)$
$\displaystyle\times\quad P_{\ell\,-m}\left(\cos
X^{\theta}_{L}\right)\,\mathfrak{T}_{m}\left(\left|X^{i}_{L}\right|,\,X^{\theta}_{L}\right)d\left|X^{i}_{L}\right|\,dX^{\theta}_{L}$

From Eq. (B), one can immediately notice that the leading order deviation for
the particular solution of the Ricci scalar for the $L^{th}$ body zone comes
at $\mathscr{O}\left(\epsilon^{4}\right)$ for $\ell\,=\,0$ with
$\mathfrak{M}_{L\,0\,0}\,\equiv\,M_{L}$ being the mass monopole, and was found
to be

$\displaystyle{}_{(4)}\mathfrak{R}_{part}$ $\displaystyle\equiv$
$\displaystyle{}_{(4)}\mathfrak{R}\,=\,-8\,\pi\,\gamma^{2}\,\epsilon^{4}\,\frac{M_{L}\,\,e^{-\frac{\gamma\,\left|Z^{i}_{L}\right|}{\sqrt{\epsilon}}}}{\left|Z^{i}_{L}\right|}$
(116)

##  Appendix C Contribution of the perturbed Ricci scalar to the dynamics of
$\mathfrak{\tilde{h}}^{\mu\nu}$, and changes to the GR two-body equations of
motion till 6 PN

The leading order particular Ricci deviation densities are given by the sum of
contributions from the two body zones

$\displaystyle{}_{(4)}\mathfrak{R}$ $\displaystyle=$
$\displaystyle-8\,\pi\,\gamma^{2}\,\epsilon^{4}\,\sum_{L=1,\,2}\,\frac{M_{L}\,\,e^{-\frac{\gamma\,\left|Z^{i}_{L}\right|}{\sqrt{\epsilon}}}}{\left|Z^{i}_{L}\right|}$
(117)

The new gothic tensor deviation trace at the leading order is given as follows

$\displaystyle{}_{(2)}\tilde{\mathfrak{h}}$ $\displaystyle\equiv$
$\displaystyle
4\,\epsilon^{2}\,\sum_{L\,=\,1,\,2}\,\frac{M_{L}}{\left|Z^{i}_{L}\right|}$
(118)

The spatial derivatives of the two scalar densities are given as follows

$\displaystyle\partial_{j}\left({}_{(2)}\mathfrak{\tilde{h}}\right)$
$\displaystyle=$
$\displaystyle-4\,\epsilon^{2}\,\sum_{L=1,\,2}\,\frac{M_{L}}{\left|Z^{i}_{L}\right|^{3}}\,Z^{L}_{j}$
(119) $\displaystyle{}_{(4)}\mathfrak{R}_{\,,\,j}$ $\displaystyle=$
$\displaystyle
8\,\pi\,\gamma^{2}\,\epsilon^{4}\,\sum_{L=1,\,2}\,\frac{M_{L}}{\left|Z^{i}_{L}\right|^{3}}\,\left(1\,+\,\frac{\gamma\left|Z^{i}_{L}\right|}{\sqrt{\epsilon}}\right)\,e^{-\frac{\gamma\,\left|Z^{i}_{L}\right|}{\sqrt{\epsilon}}}\,Z_{j}^{L}$
(120)

Therefore, the first appearance of the effective source term
$t^{\mu\nu}_{eff}$ in the $\epsilon$ series of $\tilde{\Lambda}^{\mu\nu}$,
following the definition in Eq. (82), was found to be
$\mathscr{O}\left(\epsilon^{12}\right)$ (since
$\left(-g\right)^{-2}\,\sim\,\epsilon^{4}$ and
$\tilde{\mathfrak{h}}\,\sim\,\epsilon^{2}$) in the following

$\displaystyle
t^{i\,j}_{eff}=\frac{\epsilon^{12}}{192\,\pi\,\gamma^{4}}\,\left[2\,\epsilon^{2}\left(\eta^{i\,k}\,\eta^{j\,l}\,+\,\eta^{i\,l}\,\eta^{j\,k}\,-\,\eta^{i\,j}\,\eta^{k\,l}\right)\,{}_{(4)}\mathfrak{R}\,\partial_{k}{}_{(4)}\mathfrak{R}\,\partial_{l}\left({}_{(2)}\tilde{\mathfrak{h}}\right)\right.$
$\displaystyle\left.\,+\,\left(\eta^{i\,k}\,\eta^{j\,l}\,-\,\frac{1}{2}\,\eta^{i\,j}\,\eta^{k\,l}\right)\,\left\\{4\,\partial_{k}\left({}_{(4)}\mathfrak{R}\right)\,\partial_{l}\left({}_{(4)}\mathfrak{R}\right)\,+\epsilon^{4}\,{}_{(4)}\mathfrak{R}_{hom}^{2}\,\partial_{k}\left({}_{(2)}\tilde{\mathfrak{h}}\right)\,\partial_{l}\left({}_{(2)}\tilde{\mathfrak{h}}\right)\right\\}\right]$
(121)

Eq. (118) is proportional the Newtonian potential with the ’usual definition
of mass’, and $\left|Z^{i}_{L}\right|$ are the radial distances from the
center of mass of the two compact objects to a point in the external Minkowski
space-time.

The change in the equations of motion is given by the boundary independent
part of the following surface integral, as was defined in Eq. (36)

$\displaystyle F^{i}_{eff}$ $\displaystyle=$
$\displaystyle-\epsilon^{-4}\,\oint_{\partial\mathscr{B}_{1}}\,dS_{j}\,t^{i\,j}_{eff}$
(122) $\displaystyle=$
$\displaystyle-\epsilon^{-4}\,\oint_{\partial\mathscr{B}_{1}}\,\left|dS\right|\,r^{1}_{j}\,t^{i\,j}_{eff}$
(123) $\displaystyle=$
$\displaystyle-\epsilon^{-4}\,\oint_{\partial\mathscr{B}_{1}}\,\left(\epsilon\,a_{1}\right)^{2}\,\sin\theta\,d\theta\,d\phi\,r^{1}_{j}\,t^{i\,j}_{eff}$
(124)

where $a_{1}$ is the radius of the first body zone coordinate system as seen
by an observer at asymptotic infinity, and $\epsilon\,a_{1}$ is the radius as
seen by an observer in the first body zone coordinate system. The radius is
$\epsilon\,a_{1}$ instead of just $a_{1}$ since the integral is being
performed in the body zone coordinate system to account for strong internal
gravity. $r^{1}_{k}$ is a unit one-form perpendicular at the surface of the
sphere that forms the boundary of the first body zone coordinate system.
Substituting Eq. (119) and (120) in Eq. (121), one obtains, after contracting
the extra effective energy-momentum-stress density tensor $t^{i\,j}_{eff}$
with the unit 1-form $r_{j}^{1}$

$\displaystyle
r_{j}^{i}\,t^{i\,j}_{eff}=\frac{\epsilon^{12}}{192\,\pi\,\gamma^{4}}\,\left[-512\,\pi^{2}\,\epsilon^{2}\,\gamma^{4}\,\sum_{L_{1},\,L_{2},\,L_{3}\,=\,1}^{2}\left(Z^{i}_{L_{2}}\,r_{1}^{l}\,Z_{l}^{L_{3}}\,+\,Z^{i}_{L_{3}}\,r^{l}_{1}\,Z^{L_{2}}_{l}\,-\,r^{i}_{1}\,Z^{l}_{L_{2}}\,Z^{L_{3}}_{l}\right)\,M_{L_{1}}\,M_{L_{2}}\,M_{L_{3}}\right.$
$\displaystyle\times\,\frac{\left(1\,+\,\gamma\,\left|Z^{i}_{L_{2}}\right|\right)}{\left|Z^{i}_{L_{1}}\right|\,\left|Z^{i}_{L_{2}}\right|^{3}\,\left|Z^{i}_{L_{3}}\right|^{3}}\,e^{-\frac{\gamma}{\sqrt{\epsilon}}\,\left(\left|Z^{i}_{L_{1}}\right|\,+\,\left|Z^{i}_{L_{2}}\right|\right)}\,+\,256\,\pi^{2}\,\gamma^{4}\,\sum_{L_{1},\,L_{2}\,=\,1}^{2}\,\left(r_{1}^{l}\,Z^{i}_{L_{1}}\,Z^{L_{2}}_{l}\,-\,\frac{1}{2}\,r^{i}_{1}\,Z^{l}_{L_{1}}\,Z^{L_{2}}_{l}\right)$
$\displaystyle\times\,\frac{M_{L_{1}}\,M_{L_{2}}\,\left(1+\frac{\gamma}{\sqrt{\epsilon}}\,\left|Z^{i}_{L_{1}}\right|\,\right)\,\left(1+\frac{\gamma}{\sqrt{\epsilon}}\,\left|Z^{i}_{L_{1}}\right|\,\right)}{\left|Z^{i}_{L_{1}}\right|^{3}\,\left|Z^{i}_{L_{2}}\right|^{3}}\,e^{-\frac{\gamma}{\sqrt{\epsilon}}\left(\left|Z^{i}_{L_{1}}\right|\,+\,\left|Z^{i}_{L_{2}}\right|\right)}\,+\,1024\,\pi^{2}\,\epsilon^{4}\,\gamma^{4}\,\sum_{L_{1},\,L_{2},\,L_{3}\,L_{4}\,=\,1}^{2}$
$\displaystyle\left.\times\,\left(r_{1}^{l}\,Z^{i}_{L_{3}}\,Z^{L_{4}}_{l}\,-\,\frac{1}{2}\,r^{i}_{1}\,Z^{l}_{L_{3}}\,Z^{L_{4}}_{l}\right)\,\frac{M_{L_{1}}\,M_{L_{2}}\,M_{L_{3}}\,M_{L_{4}}}{\left|Z^{i}_{L_{1}}\right|\,\left|Z^{i}_{L_{2}}\right|\,\left|Z^{i}_{L_{3}}\right|^{3}\,\left|Z^{i}_{L_{4}}\right|^{3}}\,e^{-\frac{\gamma}{\sqrt{\epsilon}}\left(\left|Z^{i}_{L_{1}}\right|\,+\,\left|Z^{i}_{L_{2}}\right|\right)}\right]$
(125)

Contribution to the surface integral of the first two summations vanish for
all combinations of $\left\\{L_{1},\,L_{2}\right\\}$ and
$\left\\{L_{1},\,L_{2},\,L_{3}\right\\}$, whereas the quartic product of the
last term has only non-vanishing contribution to the surface integral for
$\left\\{L_{1},\,L_{2},\,L_{3},\,L_{4}\right\\}\,=\,\left\\{2,\,2,\,2,\,1\right\\}\,\&\,\left\\{2,\,2,\,1,\,2\right\\}$,
which then leads to the following

$\displaystyle
F^{i}_{eff}=-\frac{32\,\pi\,\epsilon^{12}}{3}\,M_{1}\,M_{2}^{3}\,\lim\limits_{\epsilon\,\rightarrow\,0}\,\lim\limits_{a_{1}\,\rightarrow\,\infty}\oint_{\partial\mathscr{B}_{1}}\,\left(\epsilon\,a_{1}\right)^{2}\,\sin\theta\,d\theta\,d\phi\,\,\,\frac{e^{-\frac{2\,\gamma}{\sqrt{\epsilon}}\left|-r_{1\,2}\,n^{i}\,+\,\epsilon\,a_{1}\,r_{1}^{i}\right|}}{\left|-r_{1\,2}\,n^{i}\,+\,\epsilon\,a_{1}\,r_{1}^{i}\right|^{5}\,\left(\epsilon\,a_{1}\right)^{3}}$
$\displaystyle\times\,\left[r_{1}^{l}\,\left(-r_{1\,2}\,n^{i}\,+\,\epsilon\,a_{1}\,r^{i}_{1}\right)\,\epsilon\,a_{1}\,r_{l}^{1}\,-\,\frac{1}{2}\,r_{1}^{i}\,\left(-r_{1\,2}\,n^{l}\,+\,\epsilon\,a_{1}\,r^{l}_{1}\right)\,\epsilon\,a_{1}\,r_{l}^{1}\right]$
(126)
$\displaystyle=\frac{32\,\pi\,\epsilon^{12}}{3}\,\frac{M_{1}\,M_{2}^{3}}{r_{1\,2}^{4}}\,e^{-\frac{2\,\gamma\,r_{1\,2}}{\sqrt{\epsilon}}}\,\oint_{\partial\mathscr{B}_{1}}\,\left(n^{i}\,-\,\frac{1}{2}\,r^{i}_{1}\,\cos\theta\right)\,\sin\theta\,d\theta\,d\phi$
(127)
$\displaystyle=-\frac{32\,\pi\,\epsilon^{12}}{9}\,\,\frac{M_{1}\,M_{2}^{3}}{r_{1\,2}^{4}}\,e^{-\frac{2\,\gamma\,r_{1\,2}}{\sqrt{\epsilon}}}\,n^{i}$
(128)

which is an effect that comes at the 6 PN order in the equations of motion,
which is a very small deviation from GR. In the above $r_{1}^{i}$ is the unit
normal radially outward 3-vector on the surface of the sphere
$\partial\mathscr{B}_{1}$, or the body zone boundary of the first object.
$n^{i}$ is the unit vector pointing from the COM of the first body zone to the
COM of the second body zone, whereas $\theta$ is the angle between $n^{i}$ and
$r_{1}^{i}$. About $\epsilon\,\rightarrow\,0$, Eq. (C) requires $a_{1}$ or the
radius of the body zone boundary to go to infinity, as was described in [14].
Once the limit is taken, the only part of Eq. (C) that is independent of the
body zone boundary $a_{1}$, is given by Eq. (128).

The above leads to the following equations of motion till 6 PN

$\displaystyle M_{1}\,\frac{dv_{1}^{i}}{ds}$ $\displaystyle=$ $\displaystyle
F^{i}_{Newtonian}\,+\,\epsilon^{2}\,F^{i}_{1PN,GR}\,+\,\epsilon^{4}\,F^{i}_{2PN,GR}\,+\,\epsilon^{5}\,F^{i}_{2.5PN,GR}$
(129)
$\displaystyle\,+\,F^{i}_{3PN,GR}\,+\,...\,+\,\epsilon^{12}\,\left(F^{i}_{6PN,GR}\,-\,\frac{32\,\pi\,M_{1}\,M_{2}^{3}}{9\,r_{1\,2}^{4}}\,\,e^{-\frac{2\,\gamma\,r_{1\,2}}{\sqrt{\epsilon}}}\,n^{i}\right)$

where the forms of the various PN forces till
$\mathscr{O}\left(\epsilon^{6}\right)$ can be found in [52], whereas higher
order dynamics can be found in [53, *Bini:2020rzn, *Bini:2020wpo] and [56,
*Blumlein:2021txj]. Eq. (129) is an equation of motion of two classical point
particles with the mutual force given by the RHS of the above. The potential
is given by the the line integral of the the conservative acceleration/force-
per-unit-mass along a radial direction, which in this case is along $-n^{i}$
or pointing away from the COM of the first body zone (antiparallel to the
displacement vector between the first and second body zone), and hence
repulsive.

## References

  * Will [2014] C. M. Will, Living Rev. Relativity 17, 4 (2014).
  * Wildberg and Furley [2014] C. Wildberg and W. Furley,  _Philoponus: Corollaries on Place and Void with Simplicius: Against Philoponus on the Eternity of the World_, Ancient Commentators on Aristotle (Bloomsbury Publishing, 2014).
  * Galilei [2013] G. Galilei,  _Discorsi e dimostrazioni matematiche_ (Elsevier Science, 2013).
  * Newton  _et al._ [1871] I. Newton, W. Kelvin, and H. Blackburn,  _Philosophiae Naturalis Principia Mathematica_ (Innys, 1871).
  * Eötvös [1907] R. Eötvös,  _Bestimmung der Gradienten der Schwerkraft und ihrer Niveauflächen mit Hilfe der Drehwaage_ (EJ Brill, 1907).
  * Roll  _et al._ [1964] P. Roll, R. Krotkov, and R. Dicke, Annals of Physics 26, 442 (1964).
  * Schlamminger  _et al._ [2008] S. Schlamminger, K.-Y. Choi, T. A. Wagner, J. H. Gundlach, and E. G. Adelberger, Phys. Rev. Lett. 100, 041101 (2008).
  * Ransom  _et al._ [2014] S. M. Ransom, I. H. Stairs, A. M. Archibald, J. W. T. Hessels, D. L. Kaplan, M. H. van Kerkwijk, J. Boyles, A. T. Deller, S. Chatterjee, A. Schechtman-Rook, A. Berndsen, R. S. Lynch, D. R. Lorimer, C. Karako-Argaman, V. M. Kaspi, V. I. Kondratiev, M. A. McLaughlin, J. van Leeuwen, R. Rosen, M. S. E. Roberts, and K. Stovall, Nature (London) 505, 520 (2014), arXiv:1401.0535 [astro-ph.SR] .
  * Archibald  _et al._ [2018] A. M. Archibald, N. V. Gusinskaia, J. W. T. Hessels, A. T. Deller, D. L. Kaplan, D. R. Lorimer, R. S. Lynch, S. M. Ransom, and I. H. Stairs, Nature (London) 559, 73 (2018), arXiv:1807.02059 [astro-ph.HE] .
  * Damour [1982] T. Damour, in  _Les Houches Summer School on Gravitational Radiation_ (1982).
  * Damour [1984] T. Damour, Fundam. Theor. Phys. 9, 89 (1984).
  * Starobinsky [2007] A. A. Starobinsky, JETP Letters 86, 157 (2007).
  * Nojiri and Odintsov [2006] S. Nojiri and S. D. Odintsov, eConf C0602061, 06 (2006), arXiv:hep-th/0601213 .
  * Futamase and Itoh [2007] T. Futamase and Y. Itoh, Living Reviews in Relativity 10 (2007), 10.12942/lrr-2007-2.
  * LANDAU and LIFSHITZ [1975a] L. LANDAU and E. LIFSHITZ, in  _The Classical Theory of Fields (Fourth Edition)_, Course of Theoretical Physics, Vol. 2, edited by L. LANDAU and E. LIFSHITZ (Pergamon, Amsterdam, 1975) fourth edition ed., pp. 345 – 357.
  * LANDAU and LIFSHITZ [1975b] L. LANDAU and E. LIFSHITZ, in  _The Classical Theory of Fields (Fourth Edition)_, Course of Theoretical Physics, Vol. 2, edited by L. LANDAU and E. LIFSHITZ (Pergamon, Amsterdam, 1975) fourth edition ed., pp. 295 – 344.
  * LANDAU and LIFSHITZ [1975c] L. LANDAU and E. LIFSHITZ, in  _The Classical Theory of Fields (Fourth Edition)_, Course of Theoretical Physics, Vol. 2, edited by L. LANDAU and E. LIFSHITZ (Pergamon, Amsterdam, 1975) fourth edition ed., pp. 259 – 294.
  * Futamase and Schutz [1983] T. Futamase and B. F. Schutz, Phys. Rev. D 28, 2363 (1983).
  * Futamase [1987] T. Futamase, Phys. Rev. D 36, 321 (1987).
  * Itoh  _et al._ [2000] Y. Itoh, T. Futamase, and H. Asada, Phys. Rev. D 62, 064002 (2000).
  * Harte [2014] A. I. Harte, Phys. Rev. Lett. 113, 261103 (2014), arXiv:1409.4674 [gr-qc] .
  * Schwarzschild [1916] K. Schwarzschild, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys. ) 1916, 189 (1916), arXiv:physics/9905030 .
  * Fromholz  _et al._ [2014] P. Fromholz, E. Poisson, and C. M. Will, Am. J. Phys. 82, 295 (2014), arXiv:1308.0394 [gr-qc] .
  * Einstein and Infeld [1949] A. Einstein and L. Infeld, Canadian Journal of Mathematics 1, 209 (1949).
  * Infeld and Schild [1949] L. Infeld and A. Schild, Reviews of Modern Physics 21, 408 (1949).
  * Chandrasekhar [1965] S. Chandrasekhar, Astrophys. J. 142, 1488 (1965).
  * Chandrasekhar [1967] S. Chandrasekhar, Astrophys. J. 147, 334 (1967).
  * Chandrasekhar and Nutku [1969] S. Chandrasekhar and Y. Nutku, Astrophys. J. 158, 55 (1969).
  * Chandrasekhar [1969] S. Chandrasekhar, Astrophys. J. 158, 45 (1969).
  * Chandrasekhar and Esposito [1970] S. Chandrasekhar and F. P. Esposito, Astrophys. J. 160, 153 (1970).
  * Giulini [2018] D. Giulini, Int. J. Geom. Meth. Mod. Phys. 15, 1850182 (2018), arXiv:1808.09320 [math-ph] .
  * Capozziello  _et al._ [2010] S. Capozziello, A. Stabile, and A. Troisi, International Journal of Theoretical Physics 49, 1251 (2010).
  * Felice and Tsujikawa [2010] A. D. Felice and S. Tsujikawa, Living Reviews in Relativity 13 (2010), 10.12942/lrr-2010-3.
  * Nojiri and Odintsov [2011] S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59 (2011), arXiv:1011.0544 [gr-qc] .
  * Nojiri  _et al._ [2017] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * [36] H.-J. Schmidt, Astronomische Nachrichten 307, 339.
  * [37] P. Teyssandier, Astronomische Nachrichten 311, 209.
  * Cognola  _et al._ [2005] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, and S. Zerbini, JCAP 02, 010 (2005), arXiv:hep-th/0501096 .
  * Olmo [2005] G. J. Olmo, Phys. Rev. Lett. 95, 261102 (2005).
  * Berry and Gair [2011] C. P. L. Berry and J. R. Gair, Physical Review D - Particles, Fields, Gravitation and Cosmology 83 (2011), 10.1103/PhysRevD.83.104022.
  * Alves  _et al._ [2009] M. Alves, O. Miranda, and J. de Araujo, Physics Letters B 679, 401 (2009).
  * Bean  _et al._ [2007] R. Bean, D. Bernat, L. Pogosian, A. Silvestri, and M. Trodden, Phys. Rev. D 75, 064020 (2007).
  * Capozziello  _et al._ [2007] S. Capozziello, A. Stabile, and A. Troisi, Physical Review D 76 (2007), 10.1103/physrevd.76.104019.
  * Bazeia  _et al._ [2007] D. Bazeia, B. C. da Cunha, R. Menezes, and A. Petrov, Physics Letters B 649, 445 (2007).
  * Mirshekari and Will [2013] S. Mirshekari and C. M. Will, Phys. Rev. D 87, 084070 (2013), arXiv:1301.4680 [gr-qc] .
  * Kopeikin [2020] S. M. Kopeikin, Phys. Rev. D 102, 024053 (2020), arXiv:2006.08029 [gr-qc] .
  * de la Cruz-Dombriz  _et al._ [2009] A. de la Cruz-Dombriz, A. Dobado, and A. L. Maroto, Phys. Rev. D 80, 124011 (2009).
  * Nzioki  _et al._ [2014] A. M. Nzioki, R. Goswami, and P. K. S. Dunsby, (2014), arXiv:1408.0152 [gr-qc] .
  * Nashed and Nojiri [2020] G. G. L. Nashed and S. Nojiri, Phys. Rev. D 102, 124022 (2020).
  * Nashed [2018] G. G. L. Nashed, European Physical Journal Plus 133, 18 (2018).
  * Martín-García [2008] J. M. Martín-García, Computer Physics Communications 179, 597–603 (2008).
  * Itoh [2004] Y. Itoh, Phys. Rev. D 69, 064018 (2004).
  * Bini  _et al._ [2020a] D. Bini, T. Damour, and A. Geralico, Phys. Rev. D 102, 024061 (2020a), arXiv:2004.05407 [gr-qc] .
  * Bini  _et al._ [2021] D. Bini, T. Damour, A. Geralico, S. Laporta, and P. Mastrolia, Phys. Rev. D 103, 044038 (2021), arXiv:2012.12918 [gr-qc] .
  * Bini  _et al._ [2020b] D. Bini, T. Damour, and A. Geralico, Phys. Rev. D 102, 024062 (2020b), arXiv:2003.11891 [gr-qc] .
  * Blümlein  _et al._ [2021a] J. Blümlein, A. Maier, P. Marquard, and G. Schäfer, Nucl. Phys. B 965, 115352 (2021a), arXiv:2010.13672 [gr-qc] .
  * Blümlein  _et al._ [2021b] J. Blümlein, A. Maier, P. Marquard, and G. Schäfer, Phys. Lett. B 816, 136260 (2021b), arXiv:2101.08630 [gr-qc] .