# Gravitational radiation from binary systems in $f(R)$ gravity: A semi-
classical approach

Ashish Narang  Subhendra Mohanty  Soumya Jana,11footnotetext: Corresponding
Author. Email<EMAIL_ADDRESS>
###### Abstract

The rate of energy loss and orbital period decay of quasi- stable compact
binary systems are derived in $f(R)$ theory of gravity using the method of a
single vertex graviton emission process from a classical source. After
linearising the $f(R)$ action written in an equivalent scalar-tensor format in
the Einstein frame, we identify the appropriate interaction terms between the
massless spin-2 tensor mode, massive scalar mode, and the energy momentum
tensor. The definition of the scalar field is related to the $f(R)$ models.
Then using the interaction vertex we compute the rate of energy loss due to
spin-2 quadrupole radiation, which comes out to be the same as the Peter-
Mathews formula with a multiplication factor, and also the energy loss due to
the scalar dipole radiation. The total energy loss is the sum of these two
contributions. Our derivation is most general as it is applicable for both
arbitrary eccentricity of the binary orbits and arbitrary mass of the scalar
field. Using the derived theoretical formula for the period decay of the
binary systems, we compare the predictions of $f(R)$ gravity and general
relativity for the observations of three binary systems, i.e. Hulse-Taylor
Binary, PSR J1141-6545 and PSR J1738+0333. Thus we put bound on three well-
known $f(R)$ dark energy models, namely the Hu-Sawicki, the Straobinsky, and
the Tsujukawa model. We get the best constraint on $f^{\prime}(R_{0})-1$
(where $R_{0}$ is the scalar curvature of the Universe at the present epoch)
from the Tsujikawa model, i.e $|f^{\prime}(R_{0})-1|<3.44\times 10^{-4}$. This
bound is stronger than those from most of the astrophysical observations and
even some cosmological observations.

##  1 Introduction

General relativity (GR), since its inception in 1916, is the most successful
classical theory of gravity which has passed all experimental tests [1] and
more so with the recent direct detection of gravitational waves by the LIGO-
VIRGO collaboration [2] and the observation of the black hole images by the
Event Horizon Telescope [3]. However, some long-standing puzzles in GR are
still not solved and understood with full consistency. One such puzzle is
understanding the observed accelerated expansion of the Universe which
motivate many researchers to pursue modified gravity theories in the classical
framework with deviation from GR in the infrared energy scales. The simplest
modification in this direction is the $f(R)$ theory of gravity which is a
generalization of the Einstein-Hilbert action by replacing the Ricci scalar
($R$) with a function $f(R)$ (see [4, 5, 6] and the references therein for a
review). In some $f(R)$ models, the cosmological constant and the dark energy,
i.e., a new exotic form of matter are not necessary for the late time
acceleration of the Universe. Initial versions of such models [7, 8, 9, 10]
usually suffer from various instability problems [11, 12, 13, 14] and also do
not satisfy the local gravity constraints [15, 16, 17, 18]. Later, several
class of $f(R)$-dark energy models were proposed [19, 20, 21, 22, 23, 24, 25,
26] which are stable and do satisfy cosmological and Solar-System constraints
under certain limits on the parameter space. In general, $f(R)$ gravity
carries a massive scalar degree of freedom apart from the usual massless
spin-2 tensor modes [27, 28]. Dynamically, $f(R)$ gravity is equivalent to
Einstein gravity minimally coupled to a scalar field in the Einstein frame,
where the scalar field is associated with a nontrivial potential that depends
upon the form of the $f(R)$ and couples to matter through the trace of the
energy-momentum tensor. In the non-relativistic limit, the scalar field
sources a (finite-range) fifth force which is added to the usual Newtonian
force. The role of this extra scalar field in gravitational radiation and
weak-field metric for simple sources was studied in [29] using the linearized
form of $f(R)$ gravity. In some $f(R)$ theories, the fifth force can be
screened only at the galactic or Solar-System scales through the chameleon
mechanism [21, 30, 31]. This mechanism facilitates the above mentioned viable
models to conform the local gravity constraints as well as the modified
dynamics at the large scale. Recently, in Refs. [32, 33], the authors have
discussed how such screening mechanisms in scalar-tensor theories affect the
gravitational radiation from compact binary systems.

It is interesting to ask what will happen if mass of the scalar field vanishes
or more precisely, is there any vDVZ discontinuity in the $f(R)$ theory? The
quest is originally linked to the Fierz-Pauli (FP) theory of massive gravity
[34]. Generically, in a massive field theory a particle exchange gives rise to
a $(1/r)e^{-m_{g}r}$ Yukawa potential which goes to the $1/r$ potential in the
$m_{g}\rightarrow 0$ limit. The FP theory of massive graviton has the
peculiarity that in the zero graviton mass limit the Lagrangian goes smoothly
to Einstein-Hilbert (EH) linearized gravity theory, while the graviton
propagator has additional contributions from the scalar modes of the metric
which do not decouple in the zero graviton mass limit. As a result, the
Newtonian potential in the zero-mass limit of FP theory is a factor $(4/3)$
larger than the prediction from the EH theory (which of course agrees with the
Newtonian potential). The $f(R)$ theory also shares the same story. This
peculiarity of the FP theory where the action goes to the EH theory in the
zero mass limit but the graviton propagator does not was first pointed out by
van Dam and Veltman [35] and independently by Zakharov [36] and this feature
which arises in most massive gravity theories [37, 38, 39, 40] is called the
van Dam-Veltman-Zakharov (vDVZ) discontinuity (although, in the nonlinear FP
theory, a proper decoupling limit will display the vDVZ discontinuity already
in the action). However, $f(R)$ does not turn into EH action in the zero mass
limit of the scalar field unlike FP theory. Thus the anomaly of the zero mass
limit in $f(R)$ theory is different than the vDVZ discontinuity of the other
massive gravity theories. This anomaly can be avoided in some class of $f(R)$
theories which admit the chameleon screening mechanism.

Previously, gravitational radiation from compact binary systems in $f(R)$
theories has been studied [41, 42, 43]. Compact binary systems are excellent
laboratory to test theories of gravity in the strong field regime. The first
indirect evidence of Gravitational Wave (GW) radiation was obtained from
precision measurements of the Hulse -Taylor binary system [44, 45, 46]. The
orbital period loss of the compact binary system confirms Einstein’s GR [47]
to $\sim 0.1\%$ accuracy [48]. Following the Hulse-Taylor binary there have
been other precision measurements from compact binary systems [49, 50, 51].

The objective of the present paper is to derive the rate of energy loss of
binary systems (in quasi-stable orbits) due to gravitational radiation in
observationally viable $f(R)$ theories using a one graviton vertex process and
compare the predictions of GR and the predictions of $f(R)$ theories. Earlier
in [52] authors derived the rate of period decay of the binary pulsar systems
by using entirely classical method. In [53] the period decay rate was computed
classically and was used to constrain various $f(R)$ models using observed
period decay of the binary pulsar systems. However, the authors did not
consider eccentricity in their calculation. GR in the weak field limit can be
treated as a quantum field theory of spin-2 fields in the Minkowski space [54,
55, 56, 57]. Any classical gravity interaction like Newtonian potential
between massive bodies or bending of light by a massive body can be described
by a tree level graviton exchange. The result of the tree level exchange
should match the weak field classical GR results. The derivation of
gravitational radiation from binary stars as a single vertex Feynman diagram
of massless graviton emission from a classical source has been performed in
[58, 59] and the results match with the result of Peter and Mathews [47] who
used the quadrupole formula of classical GR. This method was also used in
computation of other radiations such as the vector gauge boson radiation [60]
and the massive graviton radiation [61]. In section II, the linearized $f(R)$
action is derived in the scalar-tensor form. Using the linearized action, the
rate of energy loss of a binary system due to gravitational radiation is
obtained in section III. Using the theoretical results, observational
constraints on three well known $f(R)$ dark energy models, i.e. Hu-Sawicki,
Starobinsky, and Tsujikawa models are also obtained for three NS-NS/WD binary
systems. In section V, we summarize our results.

##  2 Linearized f(R) action

The $f(R)$ gravitational action is given by,

$S=\int\sqrt{-g}\mathrm{d^{4}x}\left[-\frac{1}{16\pi
G}f(R)+\mathcal{L}_{M}\right]$ (2.1)

where $\mathcal{L}_{M}$ is the matter part of the action. Expanding the
function $f(R)$ in the Taylor series around $R=0$, we get

$f(R)=f(0)+f^{\prime}(0)R+\frac{f^{\prime\prime}(0)}{2!}R^{2}+\mathcal{O}(R^{3})+...$
(2.2)

Assuming $f(0)=0$, $f^{\prime}(0)=1$ and
$\frac{1}{2}f^{\prime\prime}(0)=\alpha$, the series becomes

$f(R)=R+\alpha R^{2}+\mathcal{O}(R^{3})+...$ (2.3)

Next we expand the gravitational action with respect to
$g_{\mu\nu}=\eta_{\mu\nu}+\kappa h_{\mu\nu}$ and $\kappa=\sqrt{32\pi G}$,
where $h_{\mu\nu}$ is the gravitation wave propagating in the flat the
background. The Ricci tensor and the Ricci scalar are expanded upto first
order in $h_{\mu\nu}$ as

$\displaystyle R_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}\left[\partial_{\mu}\partial_{\rho}h^{\rho}_{\nu}-\partial_{\nu}\partial_{\rho}h^{\rho}_{\mu}-\partial_{\mu}\partial_{\nu}h-\Box
h_{\mu\nu}\right]+\mathcal{O}(h^{2})$ (2.4) $\displaystyle R$ $\displaystyle=$
$\displaystyle
g^{\mu\nu}R_{\mu\nu}=\kappa(\partial_{\mu}\partial_{\nu}h^{\mu\nu}-\Box
h)+\mathcal{O}(h^{2}),$ (2.5)

where $\Box\equiv\partial^{\mu}\partial_{\mu}$. Therefore, the action (2.1) at
first order in $h_{\mu\nu}$ can be written as

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int\mathrm{d^{4}x}\left[\frac{1}{2}h^{\mu\nu}\Box
h_{\mu\nu}-h^{\mu\nu}\partial_{\mu}\partial^{\rho}h_{\rho\nu}-\frac{1}{2}h\Box
h+\frac{1}{2}h\partial_{\mu}\partial_{\nu}h^{\mu\nu}+\frac{1}{2}h^{\mu\nu}\partial_{\mu}\partial_{\nu}h\right.$
(2.6)
$\displaystyle+\left.2\alpha(h^{\mu\nu}\partial_{\mu}\partial_{\nu}\partial_{\alpha}\partial_{\beta}h^{\alpha\beta}+h\Box^{2}h-h^{\mu\nu}\partial_{\mu}\partial_{\nu}\Box
h-h\Box\partial_{\mu}\partial_{\nu}h^{\mu\nu})+\frac{\kappa}{2}h^{\mu\nu}T_{\mu\nu}\right]$
$\displaystyle=$
$\displaystyle\int\left[d^{4}xh_{\mu\nu}\mathcal{E}^{\mu\nu,\alpha\beta}h_{\alpha\beta}+\frac{\kappa}{2}h^{\mu\nu}T_{\mu\nu}\right],$
(2.7)

where $T^{\mu\nu}$ is the energy-momentum tensor of the source of the
gravitational waves and the kinetic operator
$\mathcal{E}^{\mu\nu,\alpha\beta}$ can be written as

$\mathcal{E}^{\mu\nu,\alpha\beta}=\left[\frac{1}{2}P^{(2)\mu\nu,\alpha\beta}-P^{(s)\mu\nu,\alpha\beta}\right]\Box+6\alpha
P^{(s)\mu\nu,\alpha\beta}\Box^{2},$ (2.8)

using the definitions of spin-2 and spin-0 projection operators given by

$\displaystyle P^{(2)}_{\mu\nu,\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\theta_{\mu\alpha}\theta_{\nu\beta}+\theta_{\mu\beta}\theta_{\nu\alpha}\right)-\frac{1}{3}\theta_{\mu\nu}\theta_{\alpha\beta},$
(2.9) $\displaystyle P^{(s)}_{\mu\nu,\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\theta_{\mu\nu}\theta_{\alpha\beta},$ (2.10)
$\displaystyle\theta_{\mu\nu}$ $\displaystyle=$
$\displaystyle\eta_{\mu\nu}-\frac{\partial_{\mu}\partial_{\nu}}{\Box}.$ (2.11)

The graviton propagator $D^{(0)}_{\mu\nu,\alpha\beta}$ is the inverse of the
kinetic operator $\mathcal{E}^{\mu\nu,\alpha\beta}$ such that

$\mathcal{E}^{\mu\nu,\alpha\beta}D^{(0)}_{\mu\nu,\alpha\beta}(x-y)=i\delta^{\mu}_{(\rho}\delta^{\nu}_{\sigma)}\delta^{4}(x-y).$
(2.12)

Without appropriate gauge condition it is not possible to inverse the kinetic
operator (Eq. (2.8)) in order to get the graviton propagator (Eq. (2.12)). To
avoid this difficulty, we take an alternative route which is described in the
following section.

###  2.1 An equivalent linearized action in Scalar-Tensor form

Using the conformal transformation between the Jordan frame metric
$g_{\mu\nu}$ and the Einstein frame metric $\tilde{g}_{\mu\nu}$:

$\displaystyle g_{\mu\nu}$ $\displaystyle=$ $\displaystyle
A^{2}(\phi)\tilde{g}_{\mu\nu},$ (2.13) $\displaystyle A^{2}(\phi)$
$\displaystyle=$ $\displaystyle\frac{1}{f^{\prime}(R)},$ (2.14) $\displaystyle
f^{\prime}(R)$ $\displaystyle=$ $\displaystyle\frac{df(R)}{dR},$ (2.15)

the $f(R)$ gravity action can be rewritten in scalar-tensor form, in the
Einstein frame,

$S=\int\sqrt{-\tilde{g}}\mathrm{d^{4}x}\left[-\frac{\tilde{R}}{16\pi
G}-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi+V(\phi)\right]+S_{M}[A^{2}(\phi)\tilde{g}_{\mu\nu},\Psi],$
(2.16)

where the scalar field $\phi$ and the potential $V(\phi)$ are identified as,

$\displaystyle\phi$ $\displaystyle=$ $\displaystyle-\sqrt{\frac{3}{16\pi
G}}\ln f^{\prime}(R),$ (2.17) $\displaystyle V(\phi)$ $\displaystyle=$
$\displaystyle\frac{Rf^{\prime}(R)-f(R)}{16\pi Gf^{\prime}(R)^{2}}.$ (2.18)

Hence the modified equation of motion for the scalar field is,

$\displaystyle\Box\phi$ $\displaystyle=$
$\displaystyle\frac{dV}{d\phi}-\frac{T}{M_{Pl}\beta(\phi)}$ (2.19)

with $\beta(\phi)=M_{Pl}\frac{d\ln A}{d\phi}$. For any $f(R)$ model
$\beta(\phi)\equiv 1/\sqrt{6}$. In the non-relativistic and static limit, the
equation (2.19) becomes

$\begin{split}\nabla^{2}\phi=&\frac{\text{d}V}{\text{d}\phi}+\frac{\rho}{\sqrt{6}M_{pl}}\\\
=&\frac{\text{d}V_{eff}}{\text{d}\phi}\end{split}$ (2.20)

where $V_{eff}(\phi;\rho)=V(\phi)+\rho\ln A(\phi)$, $\rho$ being the matter
density of the local environment of the scalar field. Therefore, the minima of
the effective potential $V_{eff}(\phi)$ depends of the mass density of the
local environment and stabilizes around those minima (denoted as
$\phi_{m}(\rho)$) accordingly. Consequently, the the mass of the scalar field
depends on the matter density of local environment such that
$m^{2}_{\phi}(\rho)=\frac{\text{d}^{2}V_{eff}}{\text{d}\phi^{2}}\Big{|}_{\phi_{m}}$.
For the dark energy models of $f(R)$ gravity, such as Hu-Sawicki, Tsujikawa,
and others, the mass of the scalar field becomes heavier in the high mass
density region and lighter in the low mass density region. Thus, the scalar
fifth force can be screened in the sufficiently dense environment. This effect
is known as the Chameleon Screening [31, 22].

Assuming a static spherically symmetric source object of constant mass density
$\rho_{0}$ and radius $r_{s}$ embedded in the homogeneous background of matter
density $\rho_{b}$, the scalar field equation can be solved by the method of
matching the interior and exterior solutions. The exterior scalar field
solution is then given by [32]

$\phi(r)\approx\phi_{\text{VEV}}-M_{Pl}\frac{GM\epsilon}{r}e^{-m_{b}r},\quad~{}r\geq
r_{s},$ (2.21)

where $\phi_{VEV}$ is the minima of the scalar field in the background far
away from the source object and $m_{b}$ is the mass of the scalar field at the
background matter density $\rho_{b}$, $M$ is the mass of the source object,
and $\epsilon$ is the screened parameter given by

$\epsilon=\frac{\phi_{\text{VEV}}-\phi_{0}}{M_{Pl}\Phi_{N}},$ (2.22)

where $\phi_{0}$ is the minima of the scalar field inside the source object
and $\Phi_{N}=GM/r_{s}$ is the Newtonian surface gravitational potential of
the source object. Usually, the matter density of the source object is much
higher as compared to the matter density of the background. In our case, we
assume compact objects embedded either in the cosmological background or in
the galactic background. In such situation, $\phi_{\text{VEV}}>>\phi_{0}$ and
therefore

$\epsilon=\frac{\phi_{\text{VEV}}}{M_{Pl}\Phi_{N}}=-\sqrt{\frac{3}{2}}\frac{\ln
f^{\prime}(R_{\text{VEV}})}{\Phi_{N}}.$ (2.23)

From Eq. (2.21) we identify the scalar charge associated to a source object as

$Q=4\pi M_{Pl}GM\epsilon.$ (2.24)

The scalar fifth force between two objects with scalar charges $Q_{a}$ and
$Q_{b}$, separated by the distance $r$ is thus given by

$\begin{split}F_{5}(r)=&-\frac{1}{4\pi}\frac{Q_{a}Q_{b}}{r^{2}}e^{-m_{b}r}\\\
=&-\frac{GM_{a}M_{b}\epsilon_{a}\epsilon_{b}}{2r^{2}}e^{-m_{b}r}.\end{split}$
(2.25)

From the Cassini mission, the solar system bound on the PPN parameter leads to
the upper bound on the screened parameter $|\epsilon|<2.3\times 10^{-5}$ [62].
For compact objects like Neutron star and White Dwarf, this upper bound on
$\epsilon$ will be more stringent as the surface gravitational potential for
compact objects is much higher than the Sun. Therefore, we neglect the effect
of fifth force (as $F_{5}/F_{N}\approx\mathcal{O}(\epsilon^{2})$) in deriving
the Keplerian orbit of the binary compact stars.

We treat the gravitational field $h_{\mu\nu}$ as the small perturbation over
the flat background such that

$\tilde{g}_{\mu\nu}=\eta_{\mu\nu}+\kappa h_{\mu\nu}.$ (2.26)

Then we linearise the action (2.16) considering one term at a time:

$\displaystyle\sqrt{-\tilde{g}}$ $\displaystyle=$ $\displaystyle
1+\frac{\kappa}{2}h+\frac{\kappa^{2}}{8}h^{2}-\frac{\kappa^{2}}{4}h^{\alpha\beta}h_{\alpha\beta}+\mathcal{O}(\kappa^{3}),$
(2.27) $\displaystyle\sqrt{-\tilde{g}}\frac{\tilde{R}}{16\pi G}$
$\displaystyle\simeq$ $\displaystyle\frac{1}{2}h^{\alpha\beta}\Box
h_{\alpha\beta}-\frac{1}{2}h\Box
h+\partial_{\alpha}h^{\alpha\beta}\partial_{\mu}h^{\mu}_{\beta}-\partial_{\alpha}h\partial_{\mu}h^{\alpha\mu}+...,$
(2.28) $\displaystyle\sqrt{-\tilde{g}}\partial_{\mu}\phi\partial^{\mu}\phi$
$\displaystyle\simeq$
$\displaystyle(1+\frac{\kappa}{2}h+...)(\eta^{\mu\nu}-\kappa
h^{\mu\nu})\partial_{\mu}\phi\partial_{\nu}\phi$ (2.29) $\displaystyle\simeq$
$\displaystyle(\partial_{\mu}\phi)^{2},$
$\displaystyle\sqrt{-\tilde{g}}V(\phi)$ $\displaystyle\simeq$
$\displaystyle\left(1+\frac{\kappa}{2}h\right)V_{\text{VEV}}+\frac{\kappa}{2\sqrt{6}}\rho_{\text{VEV}}\left(\phi-\phi_{\text{VEV}}\right)+\frac{1}{2}m^{2}_{\phi}(\phi-\phi_{\text{VEV}})^{2}+...,$

where

$\displaystyle V_{\text{VEV}}=V(\phi_{\text{VEV}}),$ (2.31)
$\displaystyle\phi_{\text{VEV}}\simeq-\sqrt{\frac{3}{16\pi G}}\ln
f^{\prime}(R_{\text{VEV}}).$ (2.32)

Note that the effective cosmological constant
$\Lambda_{eff}=V_{\text{VEV}}/M_{Pl}^{2}$,
$m^{2}_{\phi}=\frac{d^{2}V}{d\phi^{2}}\Big{|}_{\phi_{\text{VEV}}}$, and
$\bar{\phi}=\phi-\phi_{\text{VEV}}$. We have assumed flat background, instead
of the cosmological FRW background as the time scale of cosmological evolution
is much larger than the time scale of evolution of compact binary orbits.
Further, we can safely neglect the cosmological density and the cosmological
constant as compare to the astrophysical density and, therefore, the scalar
field part of the action becomes

$\displaystyle S_{\phi}$ $\displaystyle=$
$\displaystyle\int\sqrt{-\tilde{g}}\mathrm{d^{4}x}\left[\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V(\phi)\right]$
(2.33) $\displaystyle\simeq$
$\displaystyle\int\mathrm{d^{4}x}\left[\frac{1}{2}\partial_{\mu}\bar{\phi}\partial^{\mu}\bar{\phi}-\frac{1}{2}m^{2}_{\phi}\bar{\phi}^{2}\right].$

Now moving to the matter action $S_{M}[A^{2}(\phi)\tilde{g}_{\mu\nu},\Psi]$,
we have

$\displaystyle
S_{M}[A^{2}(\phi)\tilde{g}_{\mu\nu},\Psi]=\int\sqrt{-g}\mathcal{L}(g_{\mu\nu},\Psi)\mathrm{d^{4}x}.$
(2.34)

We can expand the Lagrangian in Taylor series,

$\displaystyle\mathcal{L_{M}}(g_{\mu\nu},\Psi)=\mathcal{L_{M}}(\bar{g}_{\mu\nu},\Psi)+\frac{\partial\mathcal{L_{M}}}{\partial\bar{g}_{\mu\nu}}(g_{\mu\nu}-\bar{g}_{\mu\nu})+...$
(2.35)

where $\bar{g}_{\mu\nu}$ is the background Jordon frame metric. Also expanding
$A^{2}({\phi})$ in Taylor series we can write the Jordan frame metric as,

$\displaystyle g_{\mu\nu}\simeq A^{2}_{\text{VEV}}\eta_{\mu\nu}+\kappa
A^{2}_{\text{VEV}}h_{\mu\nu}+\chi\eta_{\mu\nu}\bar{\phi}+\kappa\chi
h_{\mu\nu}\bar{\phi}+...$ (2.36)

where $\chi=\frac{d}{d\phi}A^{2}(\phi)\Big{|}_{\phi_{\text{VEV}}}$. Thus we
identify that

$\displaystyle\bar{g}_{\mu\nu}$ $\displaystyle=$ $\displaystyle
A^{2}_{\text{VEV}}\eta_{\mu\nu},$ (2.37) and, therefore,
$\displaystyle\sqrt{-g}\mathcal{L_{M}}(g_{\mu\nu},\Psi)$ $\displaystyle=$
$\displaystyle\left(A^{4}_{\text{VEV}}+\frac{d}{d\phi}A^{4}(\phi)\Big{|}_{\phi_{\text{VEV}}}\left(\phi-\phi_{\text{VEV}}\right)+...\right)\left(1+\frac{\kappa}{2}h+...\right)$
$\displaystyle\times\left[\mathcal{L}_{M}(\bar{g}_{\mu\nu},\Psi)+\frac{\partial\mathcal{L_{M}}}{\partial\bar{g}_{\mu\nu}}\left(\chi\bar{\phi}\eta_{\mu\nu}+\kappa
A^{2}_{\text{VEV}}h_{\mu\nu}+\kappa\chi\bar{\phi}h_{\mu\nu}+...\right)\right].$

Now using the definition of the energy-momentum tensor,

$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\frac{2}{\sqrt{-g}}\frac{\partial(\sqrt{-g}\mathcal{L}_{M})}{\partial
g_{\mu\nu}},$

we write the energy-momentum tensor for background Jordon frame metric,

$\displaystyle\bar{T}^{\mu\nu}=-\left[2\frac{\partial\mathcal{L}_{M}}{\partial\bar{g}_{\mu\nu}}+A^{-2}_{\text{VEV}}\eta^{\mu\nu}\mathcal{L}_{M}\right].$
(2.39)

Hence, we get,

$\displaystyle\frac{\partial\mathcal{L}_{M}}{\partial\bar{g}_{\mu\nu}}=-\frac{1}{2}\left(\bar{T}^{\mu\nu}+A^{-2}_{\text{VEV}}\eta^{\mu\nu}\mathcal{L}_{M}(\bar{g}_{\mu\nu},\Psi)\right).$
(2.40)

In the Einstein frame, the background metric is Minkowski, i.e.
$\tilde{g}_{\mu\nu}=\eta_{\mu\nu}$ and the corresponding Energy-Momentum
tensor is $\tilde{T}_{\mu\nu}$. Therefore

$\displaystyle\tilde{T}^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\left(2\frac{\partial\mathcal{L}_{M}}{\partial\eta_{\mu\nu}}+\eta^{\mu\nu}\mathcal{L}_{M}\right),$
(2.41) $\displaystyle\bar{T}^{\mu\nu}$ $\displaystyle=$ $\displaystyle
A^{-2}_{\text{VEV}}\tilde{T}^{\mu\nu},$ (2.42) and we get,
$\displaystyle\frac{\partial\mathcal{L}_{M}}{\partial\bar{g}_{\mu\nu}}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}A^{-2}_{\text{VEV}}(\tilde{T}^{\mu\nu}+\eta^{\mu\nu}\mathcal{L}_{M}).$
(2.43)

Thus we get

$\displaystyle\sqrt{-g}\mathcal{L}_{M}(g_{\mu\nu},\Psi)$ $\displaystyle\simeq$
$\displaystyle
A^{4}_{\text{VEV}}\bar{L}_{M}-\frac{1}{2}A^{2}_{\text{VEV}}\chi\tilde{T}\bar{\phi}-A^{4}_{\text{VEV}}\frac{\kappa}{2}\tilde{T}^{\mu\nu}h_{\mu\nu}+...$
(2.44)

and hence we obtain the linearised $f(R)$ action,

$\displaystyle S_{f}$ $\displaystyle=$
$\displaystyle\int\mathrm{d^{4}x}\left[-\frac{1}{2}(\partial_{\mu}h_{\nu\rho})^{2}+\frac{1}{2}(\partial_{\mu}h)^{2}-(\partial_{\mu}h)(\partial^{\nu}h^{\mu}_{\nu})+(\partial_{\mu}h_{\nu\rho})(\partial^{\nu}h^{\mu\rho})\right.$
(2.45)
$\displaystyle\left.+\frac{1}{2}A^{2}_{\text{VEV}}\chi\tilde{T}\bar{\phi}+A^{4}_{\text{VEV}}\frac{\kappa}{2}\tilde{T}^{\mu\nu}h_{\mu\nu}-\frac{1}{2}\partial_{\mu}\bar{\phi}\partial^{\mu}\bar{\phi}+\frac{1}{2}m^{2}_{\phi}\bar{\phi}^{2}\right].$

Using the definitions of $\phi$ and $\chi(\phi)$, we can rewrite the action

$\displaystyle S_{f}$ $\displaystyle=$
$\displaystyle\int\mathrm{d^{4}x}\left[-\frac{1}{2}(\partial_{\mu}h_{\nu\rho})^{2}+\frac{1}{2}(\partial_{\mu}h)^{2}-(\partial_{\mu}h)(\partial^{\nu}h^{\mu}_{\nu})+(\partial_{\mu}h_{\nu\rho})(\partial^{\nu}h^{\mu\rho})\right.$
(2.46)
$\displaystyle\left.+\frac{\kappa}{2\sqrt{6}}A^{4}_{\text{VEV}}\tilde{T}\bar{\phi}+A^{4}_{\text{VEV}}\frac{\kappa}{2}\tilde{T}^{\mu\nu}h_{\mu\nu}-\frac{1}{2}\partial_{\mu}\bar{\phi}\partial^{\mu}\bar{\phi}+\frac{1}{2}m^{2}_{\phi}\bar{\phi}^{2}\right]$

##  3 Gravitational radiation from binary system

In this section we compute the rate of energy loss of the binary system due to
gravitational radiation using the linearized action (2.46). The net
gravitational emission rate is thus sum of the massless spin-2 graviton
radiation and the scalar radiation. Therefore following the method described
in [58, 59, 61], the net emission rate is

$d\Gamma=d\Gamma^{h}+d\Gamma^{\phi}$ (3.1)

where the spin-2 graviton emission rate is given by

$\displaystyle d\Gamma^{h}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{4}A^{8}_{\text{VEV}}\sum_{\lambda=1}^{2}|T_{\mu\nu}(k^{\prime})\epsilon_{\lambda}^{\mu\nu}(k)|^{2}2\pi\delta(\omega-\omega^{\prime})\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega}$
(3.2) $\displaystyle=$
$\displaystyle\frac{\kappa^{2}A^{8}_{\text{VEV}}}{8(2\pi)^{2}}\sum_{\lambda=1}^{2}(T_{\mu\nu}(k^{\prime})T^{*}_{\alpha\beta}(k^{\prime})\epsilon_{\lambda}^{\mu\nu}(k)\epsilon^{*\alpha\beta}_{\lambda}(k))\frac{d^{3}k}{\omega}\delta(\omega-\omega^{\prime}).$

Here $T_{\mu\nu}(k^{\prime})$ is the stress-energy tensor in the momentum
space of the spin-2 gravitons. Using the polarization sum of massless spin-2
gravitons

$\displaystyle\sum_{\lambda=1}^{2}\epsilon_{\mu\nu}^{\lambda}(k)\epsilon_{\alpha\beta}^{*\lambda}(k)=\frac{1}{2}(\eta_{\mu\alpha}\eta_{\nu\beta}+\eta_{\mu\alpha}\eta_{\nu\alpha})-\frac{1}{2}\eta_{\mu\nu}\eta_{\alpha\beta},$
(3.3)

we get

$\displaystyle
d\Gamma^{h}=\frac{\kappa^{2}A^{8}_{\text{VEV}}}{(2\pi)^{2}}\frac{\pi}{5}\left(T_{ij}(\omega^{\prime})T^{*}_{ji}(\omega^{\prime})-\frac{1}{3}|T^{i}_{i}(\omega^{\prime})|^{2}\right)\delta(\omega-\omega^{\prime})\omega\delta\omega.$
(3.4)

Therefore , the rate of energy loss due to spin-2 graviton radiation is

$\displaystyle\frac{dE^{h}}{dt}=\frac{\kappa^{2}A^{8}_{\text{VEV}}}{20\pi}\int\left(T_{ij}(\omega^{\prime})T^{*}_{ji}(\omega^{\prime})-\frac{1}{3}|T_{i}^{i}(\omega^{\prime})|^{2}\right)\omega^{2}\delta(\omega-\omega^{\prime})d\omega.$
(3.5)

The classical energy-momentum tensor for the binary system orbiting in the
$x-y$ plane is [58, 59, 61]

$\displaystyle T_{\mu\nu}(x^{\prime})$ $\displaystyle=$
$\displaystyle\mu\delta^{3}(\vec{x^{\prime}}-\vec{x}(t))U_{\mu}U_{\nu},$ (3.6)
$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\frac{M_{1}M_{2}}{M_{1}+M_{2}},\ \text{and}\
U_{\mu}=(1,\dot{x},\dot{y},0),$ (3.7)

where $M_{1}$ and $M_{2}$ are the masses of the binary stars, $\mu$ is the
reduced mass, $U_{\mu}$ is the non-relativistic four velocity of the reduced
mass.

We can write the Keplerian orbit in the parametric form as

$x=a(\cos\xi-e),\hskip 11.38092pty=a\sqrt{(1-e^{2})}\sin\xi,\hskip
11.38092pt\Omega t=\xi-e\sin\xi,$ (3.8)

where $a$ and $e$ are the semi-major axis and eccentricity of the elliptic
orbit respectively. Since the angular velocity of an eccentric orbit is not
constant, we can write the Fourier transform of the current density in terms
of the $n$ harmonics of the fundamental frequency
$\Omega=\Big{[}G\frac{(M_{1}+M_{2})}{a^{3}}\Big{]}^{\frac{1}{2}}$. Using Eq.
(3.8) we can write the Fourier transforms of the velocity components in the
Kepler orbit as

$\dot{x}_{n}=\frac{1}{T}\int^{T}_{0}e^{in\Omega t}\dot{x}dt=-ia\Omega
J^{\prime}_{n}(ne),$ (3.9)

and

$\dot{y}_{n}=\frac{1}{T}\int^{T}_{0}e^{in\Omega
t}\dot{y}dt=\frac{a\sqrt{(1-e^{2})}}{e}\Omega J_{n}(ne),$ (3.10)

where we have used $T=2\pi/\Omega$ and the Bessel function identity
$J_{n}(z)=\frac{1}{2\pi}\int^{2\pi}_{0}e^{i(n\xi-z\sin\xi)}d\xi$. The prime
over the Bessel function denotes the derivative with respect to the argument.
Hence the Fourier transforms of the orbital coordinates become

$x_{n}=\frac{\dot{x}_{n}}{-i\Omega n}=\frac{a}{n}J^{\prime}_{n}(ne),\hskip
11.38092pty_{n}=\frac{\dot{y}_{n}}{-i\Omega
n}=\frac{ia\sqrt{1-e^{2}}}{ne}J_{n}(ne).$ (3.11)

The Fourier transforms of components of the stress energy tensor with
$\omega^{\prime}=n\Omega$ are given by [59, 61]

$\displaystyle
T_{ij}(\omega^{\prime})T^{ij*}(\omega^{\prime})-\frac{1}{3}|T^{i}_{i}(\omega^{\prime})|^{2}=4\mu^{2}\omega^{4}a^{4}f(n,e),$
(3.12)

where

$\begin{split}f(n,e)=\frac{1}{32n^{2}}\Big{\\{}[J_{n-2}(ne)-2eJ_{n-1}(ne)+2eJ_{n+1}(ne)+\frac{2}{n}J_{n}(ne)-J_{n+2}(ne)]^{2}+\\\
(1-e^{2})[J_{n-2}(ne)-2J_{n}(ne)+J_{n+2}(ne)]^{2}+\frac{4}{3n^{2}}J^{2}_{n}(ne)\Big{\\}}.\end{split}$
(3.13)

Therefore,

$\displaystyle\frac{dE^{h}}{dt}$ $\displaystyle=$
$\displaystyle\frac{32G}{5}A^{8}_{\text{VEV}}\sum_{n=1}^{\infty}(n\Omega)^{2}\mu^{2}a^{4}(n\Omega)^{4}f(n,e)$
(3.14) $\displaystyle=$
$\displaystyle\frac{32G}{5}A^{8}_{\text{VEV}}\mu^{2}a^{4}\Omega^{6}(1-e^{2})^{-7/2}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right).$

Thus the emission rate of spin-2 gravitons is exactly same as the Peter-
Mathews formula [47] with the multiplication factor $A^{8}_{\text{VEV}}$. For
the scalar part of the gravitational radiation we identify the scalar
interaction Lagrangian in the linearized action (2.46). In the non-
relativistic limit the trace of the energy-momentum tensor $\tilde{T}=-\rho$,
where $\rho$ is the energy density of the binary systems. Since $\rho$
contains the scalar charges of the source objects they interact with the
scalar field and radiate energy. Therefore, the effective interaction
Lagrangian is

$\begin{split}\mathcal{L}_{s}=&\frac{\kappa
A^{4}_{\text{VEV}}}{2\sqrt{6}}\bar{\phi}\rho\\\ \simeq&\frac{\kappa
A^{4}_{\text{VEV}}}{2}\bar{\phi}\rho_{s}M_{Pl},\end{split}$ (3.15)

where the scalar charge density $\rho_{s}$ is related to the energy density
$\rho$ as $\rho\equiv\sqrt{6}M_{Pl}\rho_{s}$. The scalar charge density
$\rho_{s}(x)$ for the binary stars (denoted by $a=1,\,2$) may be written as

$\rho_{s}(x)=\sum_{a=1,2}Q_{a}\delta^{3}(\vec{x}-\vec{x_{a}}(t)),$ (3.16)

where $Q_{a}$ is the scalar charge [Eq. (2.24)] in the star and
$\vec{x_{a}}(t)$ represents the Keplerian orbit of the binary stars. Thus the
scalar part of the gravitational radiation is

$\displaystyle d\Gamma^{\phi}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}A^{8}_{\text{VEV}}}{4}M_{Pl}^{2}\int|\rho_{s}(\omega^{\prime})|^{2}(2\pi)\delta(\omega-\omega^{\prime})\frac{d^{3}k}{(2\pi)^{3}(2\omega)},$
(3.17)

where $d^{3}k=k^{2}dkd\Omega_{k}$. Using the dispersion relation for the
scalar field is $k^{2}=\omega^{2}-m^{2}_{\phi}$, we get

$d\Gamma^{\phi}=\frac{\kappa^{2}A^{8}_{\text{VEV}}M_{Pl}^{2}}{32\pi^{2}}\int|\rho_{s}(\omega^{\prime})|^{2}\omega\delta(\omega-\omega^{\prime})\sqrt{1-m_{\phi}^{2}/{\omega^{2}}}\,d\omega\,d\Omega_{k},$
(3.18)

where $\rho_{s}(\omega^{\prime})$ is the Fourier transform of the charge
density $\rho_{s}(s)$ (Eq. (3.16)) given by

$\rho_{s}(\omega^{\prime})=\int\frac{1}{T}\int^{T}_{0}e^{i\vec{k}\cdot\vec{x}}e^{-i\omega^{\prime}t}\sum_{a=1,2}Q_{a}\delta^{3}(\vec{x}-\vec{x}_{a}(t))d^{3}xdt,$
(3.19)

where $T=2\pi/\Omega$. Using the Taylor’s series expansion
$e^{i\vec{k}\cdot\vec{x}}=1+i\vec{k}\cdot\vec{x}+....$ and keeping only the
leading order contribution, we get

$\rho_{s}(\omega^{\prime})=(Q_{1}+Q_{2})\delta(\omega^{\prime})+i\mu\left(\frac{Q_{1}}{M_{1}}-\frac{Q_{2}}{M_{2}}\right)\left(k_{x}x(\omega^{\prime})+k_{y}y(\omega^{\prime})\right)+\mathcal{O}(k^{2}).$
(3.20)

The first term does not contribute in the radiation formula. Using the Fourier
transform $x(\omega^{\prime})$ and $y(\omega^{\prime})$ [Eq. (3.11)] and the
angular average
$<k_{x}^{2}>=<k_{y}^{2}>=\frac{1}{3}(n\Omega)^{2}(1-n_{0}^{2}/n^{2})$, where
$n_{0}=m_{\phi}/\Omega$, we obtain

$|\rho_{s}(\omega)|^{2}=\frac{\mu^{2}}{3}\left(\frac{Q_{1}}{M_{1}}-\frac{Q_{2}}{M_{2}}\right)^{2}\Omega^{2}a^{2}\left(1-\frac{n_{0}^{2}}{n^{2}}\right)\left[J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}^{2}(ne)\right].$
(3.21)

Using Eq. (3.21) in Eq. (3.18) we obtain the rate of change of energy loss due
to the scalar dipole radiation

$\begin{split}\frac{dE_{s}}{dt}=&\frac{\kappa^{2}M_{Pl}^{2}A^{8}_{\text{VEV}}}{24\pi}\left(\frac{Q_{1}}{M_{1}}-\frac{Q_{2}}{M_{2}}\right)^{2}\mu^{2}\Omega^{4}a^{2}\sum_{n=1}^{\infty}n^{2}\left(1-\frac{n_{0}^{2}}{n^{2}}\right)^{3/2}\left[J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}^{2}(ne)\right]\\\
=&\frac{GA^{8}_{\text{VEV}}}{3}(\epsilon_{1}-\epsilon_{2})^{2}\mu^{2}\Omega^{4}a^{2}\sum_{n=1}^{\infty}n^{2}\left(1-\frac{n_{0}^{2}}{n^{2}}\right)^{3/2}\left[J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}^{2}(ne)\right],\end{split}$
(3.22)

where we used the definition of the scalar charge $Q$ as mentioned in Eq.
(2.24). Further, using the definition of the screened parameter $\epsilon$, we
rewrite the Eq. (3.22) as

$\begin{split}\frac{dE_{s}}{dt}=&\frac{GA^{8}_{\text{VEV}}}{2}\left(\frac{\ln
f^{\prime}(R_{\text{VEV}})}{\Phi_{N,2}}-\frac{\ln
f^{\prime}(R_{\text{VEV}})}{\Phi_{N,1}}\right)^{2}\mu^{2}\Omega^{4}a^{2}\\\
&\times\sum_{n=1}^{\infty}n^{2}\left(1-\frac{n_{0}^{2}}{n^{2}}\right)^{3/2}\left[J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}^{2}(ne)\right].\end{split}$
(3.23)

For dark energy models, the mass of the scalar field $m_{\phi}\sim 10^{-33}$
eV at the cosmological scales and for the compact binary systems the typical
value of $\Omega\sim 10^{-19}$ eV. Therefore, if $m_{\phi}<<\Omega$, i.e.
$n_{0}<<1$ and $\left(1-\frac{n_{0}^{2}}{n^{2}}\right)^{3/2}\sim 1$, then
using the identity
$\sum_{n=1}^{\infty}n^{2}\left[J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}^{2}(ne)\right]=(1/4)(2+e^{2})(1-e^{2})^{-5/2}$
[47], we can rewrite the Eq. (3.23) in a more compact form

$\frac{dE_{s}}{dt}=\frac{GA^{8}_{\text{VEV}}}{4}\left(\frac{\ln
f^{\prime}(R_{\text{VEV}})}{\Phi_{N,2}}-\frac{\ln
f^{\prime}(R_{\text{VEV}})}{\Phi_{N,1}}\right)^{2}\mu^{2}\Omega^{4}a^{2}\left[\frac{1+e^{2}/2}{(1-e^{2})^{5/2}}\right].$
(3.24)

However, we use Eq. (3.23) to analyse the observation and constraining the
theory as the dipole radiation can also occur even when $m_{\phi}<\Omega$ but
$m_{\phi}>>10^{-33}$ eV at the much smaller astrophysical length scale (such
as the galactic scale). The orbital period decay is given by,

$\dot{P}_{b}=-6\pi
G^{-3/2}(M_{1}M_{2})^{-1}(M_{1}+M_{2})^{-1/2}a^{5/2}\dot{E}$ (3.25)

where $\dot{E}=\dot{E}_{h}+\dot{E}_{s}$.

##  4 Constraints from observations

In this section, using the formulae derived in the previous section we shall
compare the observational period decay of the compact binary systems (neutron
star- neutron star/white dwarf) with that predicted theoretically in three
$f(R)$ gravity models, namely, the Hu-Sawicki model [22], the Starobinsky [23]
model and the Tsujikawa model [24]. We constrain the parameters of these three
models and compare these constraints with the solar system constrains [63]. We
use three binary objects for our analysis, i.e., PSR B1913+16 (Hulse-Taylor
binary) [48], PSR J1141-6545 (High-eccentric NS-WD binary) [64, 65, 66] and
PSR J1738+0333 (Low-eccentric NS-WD binary) [51] listed in Table 1.

We define relative ‘change’ of intrinsic (i.e. observed) orbital period decay
and $f(R)$ predicted orbital period decay with respect to that predicted from
GR as,

$\displaystyle\Delta_{\text{Obs}}$ $\displaystyle=$
$\displaystyle\left|\frac{\dot{P}_{b,\text{Intrinsic}}-\dot{P}_{b,\text{GR}}}{\dot{P}_{b,\text{GR}}}\right|,$
(4.1) $\displaystyle\Delta_{\text{f(R)}}$ $\displaystyle=$
$\displaystyle\left|\frac{\dot{P}_{b,\text{f(R)}}-\dot{P}_{b,\text{GR}}}{\dot{P}_{b,\text{GR}}}\right|.$
(4.2)

The allowed $f(R)$ models must satisfy the condition
$\frac{\Delta_{f(R)}}{\Delta_{Obs}}<1$.

Parameters | PSR B1913+16 | PSR J1738+0333 | PSR J1141-6545
---|---|---|---
Pulsar mass $m_{1}$ ($M_{\odot}$) | $1.438\pm 0.001$ | $1.46^{+0.06}_{-0.05}$ | $1.27\pm 0.01$
Companion mass $m_{2}$ ($M_{\odot}$) | $1.390\pm 0.001$ | $0.181^{+0.008}_{-0.007}$ | $1.02\pm 0.01$
Eccentricity $e$ | $0.6171340(4)$ | $(3.4\pm 1.1)\times 10^{-7}$ | $0.171884(2)$
Orbital period $P_{b}$ (d) | $0.322997448918(3)$ | $0.3547907398724(13)$ | $0.1976509593(1)$
Intrinsic $\dot{P_{b}}(10^{-12}\rm{ss^{-1}})$ | $-2.398\pm 0.004$ | $(-25.9\pm 3.2)\times 10^{-3}$ | $-0.403(25)$
GR $\dot{P_{b}}(10^{-12}\rm{ss^{-1}})$ | $-2.40263\pm 0.00005$ | $-27.7^{+1.5}_{-1.9}\times 10^{-3}$ | $-0.386$
Table 1: Summary of the measured orbital parameters and the orbital period
derivative values from observation and GR for PSR B1913+16 [48], PSR
J1738+0333 [51] and PSR J1141-6545 [65]. The uncertainties in the last digits
are quoted in the parenthesis.

###  4.1 Hu-Sawicki Model

The Hu-Sawicki dark-energy model [22] is given by

$f(R)=R-m^{2}\frac{c_{1}\left(R/m^{2}\right)^{n}}{c_{2}\left(R/m^{2}\right)^{n}+1},$
(4.3)

where $n>0$, $c_{1}>0$, and $c_{2}>0$ are dimensionless parameters and $m^{2}$
is the mass scale in the theory defined by $m^{2}=8\pi G\bar{\rho}_{0}/3$,
$\bar{\rho}_{0}$ being the present day average matter density of the Universe.
The theory is designed such that it effectively approaches the $\Lambda CDM$
model in the high curvature/ redshift limit and can also explain the present
day observed accelerated expansion of the Universe, i.e.
$\lim_{R\to\infty}f(R)$ $=R-2\Lambda_{eff}$ and $\lim_{R\to 0}f(R)=0$ . The
second condition is satisfied if $n>0$. For large $R$, the model (4.3) can be
approximated as

$f(R)=R-m^{2}\frac{c_{1}}{c_{2}}+m^{2}\frac{c_{1}}{c_{2}^{2}}\left(\frac{m^{2}}{R}\right)^{n}+\mathcal{O}(R^{-2n}).$
(4.4)

Then the effective cosmological constant is identified as
$\Lambda_{eff}=\frac{m^{2}c_{1}}{2c_{2}}$ and

$\frac{c_{1}}{c_{2}}=6\frac{1-\Omega_{m0}}{\Omega_{m0}},$ (4.5)

where $\Omega_{m0}=\bar{\rho}_{0}/\rho_{c0}$ is the ratio of the present day
matter density and the present day total energy density of the Universe.

With the help of the constraint equation (4.5), we have now two free
parameters in the theory. We choose $f^{\prime}(R_{\text{VEV}})$ (i.e. the
derivative of $f(R)$ with respect to $R$ at the background of the binary
system) and $n$. Assuming the local background of the binary system is the
galactic medium where the effective scalar potential $V_{eff}(\phi)$ has a
minimum, we have $R_{\text{VEV}}\equiv R_{g}=8\pi G\rho_{g}$, where
$\rho_{g}=10^{-24}$ $gm$ $cm^{-3}$ is the mean galactic density. The mass of
the scalar field at the galactic background is given by

$m^{2}_{\phi}(\rho_{g})=\frac{1}{3}\left[\frac{1}{f^{\prime\prime}(R_{g})}-\frac{4f(R_{g})}{f^{\prime}(R_{g})^{2}}+\frac{R_{g}}{f^{\prime}(R_{g})}\right],$
(4.6)

where

$\displaystyle f(R_{g})$ $\displaystyle=$ $\displaystyle
R_{g}\left[1-2\left(\frac{1}{\Omega_{m0}}-1\right)\frac{\bar{\rho}_{0}}{\rho_{g}}-\frac{f^{\prime}(R_{g})-1}{n}\right],$
(4.7) $\displaystyle f^{\prime\prime}(R_{g})$ $\displaystyle=$
$\displaystyle-\frac{n+1}{R_{g}}\left[f^{\prime}(R_{g})-1\right].$ (4.8)

Eqs. (4.6), (4.7), and (4.8) are obtained from Eqs.(2.17), (2.18), and (4.4).
Thus the formulae for rate of energy loss due to scalar dipole radiation
(3.23), the orbital period decay (3.25), and the relative change in period
decay (4.2) can be written as the function of the free parameters
$f^{\prime}(R_{g})$ and $n$. We constrain the the parameter space
($f^{\prime}(R_{g}),n$) from the observations by imposing the conditions
$0<m^{2}_{\phi}(\rho_{g})<\Omega^{2}$ and $\Delta_{f(R)}<\Delta_{\text{Obs}}$.
In Fig. 1, we show the constraint on the parameter space. We see that
constraint on the Hulse-Taylor binary pulsar ($|f^{\prime}(R_{g})-1|<10^{-3}$)
is the weakest among all three systems we considered, where as for other two
NS-WD systems the constraints are stronger. For the high eccentric NS-WD
binary PSR J1141-6545, $|f^{\prime}(R_{g})-1|<10^{-6}$ and for the low
eccentric NS-WD binary PSR J1738+0333, $|f^{\prime}(R_{g})-1|<10^{-7}$.
However all the constraints are weaker than the Solar-System bound
$|f^{\prime}(R_{g})-1|<10^{-10}$ which comes from the Cassini Mission. We note
that all the bounds are very weakly dependent on $n$. The reason is that $n$
appears in the dipole radiation formula only through the definition of
$m_{\phi}$, where as the dominating factor is the difference in the charge to
mass ratio which does not depend on $n$.

Figure 1: ($a$): Constraints on galactic scale amplitude $f^{\prime}(R_{g})$
is presented in the parameter space $f^{\prime}(R_{g})-1$ (denoted as
$f_{R_{g}}-1$) vs. $n$ for the Hu- Sawicki model. The shaded regions are
allowed by the corresponding binary stars system. The black dashed line
indicates the boundary below which any parameter values are allowed by the
Solar-System [22]. ($b$): Constrains on cosmological scale amplitude
$f^{\prime}(R_{0})$ are shown in the $(f^{\prime}(R_{0})-1)$ vs. $n$ plane. In
the figure, $f^{\prime}(R_{0})$ is labeled as $f_{R_{0}}$.

We map the constraint parameter space ($f^{\prime}(R_{g}),n$) at the galactic
scale to the parameter space ($(f^{\prime}(R_{0})),n$), where $R_{0}$ is the
curvature of the cosmological background at the present epoch. The mapping is
done using the relation

$\frac{f^{\prime}(R_{g})-1}{f^{\prime}(R_{0})-1}=\left(\frac{R_{g}}{R_{0}}\right)^{n+1}=\left[\frac{\rho_{g}}{\bar{\rho}_{0}(4/\Omega_{m0}-3)}\right]^{n+1}.$
(4.9)

In Fig. 1, we show the constraint on the parameters translated onto the
cosmological scale. The shaded regions are allowed ones for each of the binary
systems. The lower boundary corresponds to the fact that $m_{\phi}=\Omega$.
The Solar-System constraint is bounded below the black dashed line. Here, we
note $n$ dependence in the figure.

###  4.2 Starobinsky dark energy model

In the same year of publication of Hu-Sawicki model, Starobinsky also proposed
a $f(R)$ dark energy model [23]

$f(R)=R+\lambda
R_{c}\left[\left(1+\frac{R^{2}}{R_{c}^{2}}\right)^{-k}-1\right],$ (4.10)

where $k,\lambda>0$ and $R_{c}$ is the curvature scale of the order of the
cosmological constant. For large $R$, the model can be approximated as

$f(R)=R-\lambda R_{c}+\lambda
R_{c}\left(\frac{R_{c}}{R}\right)^{2k}+\mathcal{O}(R^{-2k-2}).$ (4.11)

Therefore, we identify the effective cosmological constant in the model as
$\Lambda_{eff}=\lambda R_{c}/2$. The approximated forms of the Hu-Sawicki
model (4.4) and the Starobinsky model (4.11) are the same with the one to one
mapping:

$2k\to n,\quad~{}R_{c}\to C_{2}^{-1/n},\quad~{}\lambda\to
m^{2}C_{1}C_{2}^{1/n-1}.$

However, this mapping does not hold for full theory and is only relevant at
high curvature regime. For our purpose it is sufficient to assume the large
$R$ approximated forms and therefore we get exactly same constraint on the
Starobinsky model as we get in the Hu-Sawicki model. The only difference is
that we need to keep the mapping in the mind.

###  4.3 Tsujikawa model

In the next year of publications of the Hu-Sawicki and the Starobinsky model,
another well-known $f(R)$ dark energy model was proposed by Tsujikawa [24].
The model is given by

$f(R)=R-\lambda R_{c}\tanh{\frac{R}{R_{c}}},$ (4.12)

where $\lambda>0$ and $R_{c}>0$ are only two parameters in the theory. In the
limit $R\to\infty$, this theory approaches $\lim_{R\to\infty}f(R)=R-\lambda
R_{c}$. Therefore, here also the effective cosmological constant becomes
$\Lambda_{eff}=\lambda R_{c}/2$. With this constraint relation, the theory has
only one free parameter. We choose $\lambda$ to be the free parameter and
$R_{c}=2\Lambda_{eff}/\lambda$.

The mass of the scalar field at the galactic scale is given by Eq. (4.6),
where $f(R_{g})$, $f^{\prime}(R_{g})$, and $f^{\prime\prime}(R_{g})$ are now
given by

$\displaystyle f(R_{g})$ $\displaystyle=$ $\displaystyle
R_{g}\left[1-\frac{2\bar{\rho}_{0}}{\rho_{g}}\left(\frac{1}{\Omega_{m0}}-1\right)\tanh{\left[\frac{\lambda\rho_{g}}{2\bar{\rho}_{0}(1/\Omega_{m0}-1)}\right]}\right],$
(4.13) $\displaystyle f^{\prime}(R_{g})$ $\displaystyle=$ $\displaystyle
1-\lambda\,\text{sech}^{2}\left[\frac{\lambda\rho_{g}}{2\bar{\rho}_{0}(1/\Omega_{m0}-1)}\right],$
(4.14) $\displaystyle f^{\prime\prime}(R_{g})$ $\displaystyle=$
$\displaystyle-\frac{\lambda\rho_{g}}{R_{g}\bar{\rho}_{0}(1/\Omega_{m0}-1)}\tanh{\left[\frac{\lambda\rho_{g}}{2\bar{\rho}_{0}(1/\Omega_{m0}-1)}\right]}\left[f^{\prime}(R_{g})-1\right].$
(4.15)

By imposing the conditions $0<m^{2}_{\phi}(\rho_{g})<\Omega^{2}$ and
$\Delta_{f(R)}<\Delta_{\text{Obs}}$ we get the constraint on the parameter
$\lambda$ which can be translated into the constraint on $f^{\prime}(R_{g})-1$
, using the equation (4.14). Further the constraint can be translated to the
cosmological scale by using the relation

$f^{\prime}(R_{0})=1-\lambda\,\text{sech}^{2}\left[\frac{\lambda(4/\Omega_{m0}-3)}{2(1/\Omega_{m0}-1)}\right].$
(4.16)

In the Table 2, we show the constraints on the parameters for three binary
systems. We note that for the binary systems the constraints are of the same
order of magnitude. We get $10^{-16}<|f^{\prime}(R_{g})-1|<10^{-6}$ at the
galactic scale and $10^{-11}<|f^{\prime}(R_{0})-1|<10^{-4}$ at the
cosmological scale.

System |  $\lambda$ |  $|f_{R_{g}}-1|$ | $|f_{R_{0}}-1|$
---|---|---|---
Solar-System |  $\lambda<4.9\times 10^{-11}$ or $\lambda>1.11\times 10^{-4}$ |  $<4.9\times 10^{-11}$ | $<0.2$
Hulse-Taylor Binary |  [$1.02\times 10^{-11}$, $2.09\times 10^{-4}$] |  $[6.74\times 10^{-16},6.21\times 10^{-6}]$ | $[1.02\times 10^{-11},2.09\times 10^{-4}]$
PSR J1141-6545 |  [$7.41\times 10^{-12}$, $2.09\times 10^{-4}$] |  $[6.74\times 10^{-16},3.39\times 10^{-6}]$ | $[7.41\times 10^{-12},2.09\times 10^{-4}]$
PSR J1738+0333 |  [$1.10\times 10^{-11}$, $2.09\times 10^{-4}$] |  $[6.74\times 10^{-16},8.71\times 10^{-6}]$ | $[1.10\times 10^{-11},2.09\times 10^{-4}]$
Table 2: Bounds on the Tsujikawa Model parameter for the three binary systems.
The Bound on the Solar-System comes from the Cassini Mission [22].

##  5 Conclusions

In this paper we derive the most general formula for the rate of energy loss
of quasi-stable compact binary systems in $f(R)$ theories of gravity using a
single vertex graviton emission process from a classical source. Then we use
it to put constraints on three well known $f(R)$ dark energy models, namely
the Hu-Sawicki, the Starobinsky, and the Tsujikawa models from binary pulsar
observations.

In $f(R)$ gravity, an extra massive scalar mode appears apart from the
massless spin-2 modes. This extra scalar mode affects the orbital motion of
the binary stars in two ways. One is that an attractive short ranged “fifth"
force adds up to the usual Newtonian gravitational force between two compact
objects. The other effect is that the scalar dipole radiation carries away
some part of the total mechanical energy of the binary system. The scalar
fifth force depends in the product of the scalar charges of the binary stars
and this eventually leads to the ratio of the scalar fifth force to the
Newtonian force being proportional to the product of the screened parameters
($\epsilon_{a}\epsilon_{b}$) for the binary stars. The screened parameter is a
$f(R)$-model dependent parameter and bounded by $|\epsilon|<10^{-5}$ for the
Solar-System. The bound becomes smaller and smaller for the compact objects
like white-dwarf and neutron stars as the surface gravitational potential for
these objects are much higher than that of the Sun. Therefore, we can safely
assume that the quasi-stable orbits of binary compact stars are unaffected by
the scalar fifth force. On the other hand, the scalar dipole radiation depends
on the square of the difference of the scalar charge to mass ratio of the
binary stars, which means effectively the square of the difference in the
screened parameters (i.e. $(\epsilon_{a}-\epsilon_{b})^{2}$ ). For the binary
neutron star systems this effect could be negligible although non-vanishing if
there is a slight difference in masses of the companion stars. However, for
the neutron star- white dwarf or neutron star - black hole systems there could
be significant effect of the scalar dipole radiation. Further it is to be
noted that contribution of the scalar dipole radiation rapidly falls off when
the scalar mass is greater than the orbital period of the binary system ( i.e.
$m_{\phi}>\Omega$). Again the scalar mass depends on the $f(R)$ model and also
the background matter density where the effective scalar potential attains a
minimum.

Considering all these, we have derived the rate of energy loss due to scalar
dipole radiation from the linearized $f(R)$ gravity when written in the
scalar-tensor format in the Einstein frame. We have obtained the most general
formula, i.e. Eq. (3.23), which applies to arbitrary eccentricity of the
binary orbit and arbitrary scalar mass. This formula is different from earlier
results in the literature [52, 53], which are applicable under specific
limits. Moreover our method is semi-classical while theirs are fully
classical.

Using the general formula for the rate of energy loss, we compare the
prediction for the period decay of the binary orbits with that of the Peter-
Mathews formula using GR. We have applied our theoretical result to constrain
$f(R)$ dark energy models from the observations of the period decay of three
binary systems of varied eccentricity and types, which are PSR B1913+16
(Hulse-Taylor), PSR J1141-6545 (high eccentric NS-WD), and PSR J1738+0333 (low
eccentric NS-WD). To summarize our results, at the galactic mean density as
the background of the binary system, we found that for Hu-Sawicki model the
NS-WD systems give stronger constraints on the Hulse-Taylor binary pulsar
which is a NS-NS system. However all these contraints are much weaker than the
Solar-System constraint. Although at the Galactic scale the constraints on
$f^{\prime}(R_{g})-1$ weakly depend on the model parameter $n$, at the
cosmolgical scale the translated constraint on $f^{\prime}(R_{0})-1$ strongly
depends on $n$. The Starobinsky model has the one to one mapping with the Hu-
Sawicki model at the high curvature regime which include our cases in this
paper. The Tsujikawa model provides better constraints both at the galactic
scale and when translated onto the cosmological scale. One important advantage
of the Tsujikawa model over other two is that effectively it has only one free
parameter when the theory goes to the $\Lambda CDM$ limit in the high
curvature regime and therefore, this model is relatively more robust. Also,
from comparison with the observations, it has been noted that all three binary
systems provide similar constraints for the Tsujikawa model. The constraints
are $|f^{\prime}(R_{g})-1|<10^{-6}$ at the galactic scale and
$|f^{\prime}(R_{0})-1|<10^{-4}$ when translated at the cosmological scale. In
Table 3, we compare our result with other constraints from different
observations at different scales. We note that our bound is stronger than
those coming from most of the astrophysical observations and even some
cosmological observations, like the CMB spectrum.

Observations |  $|f^{\prime}(R_{0})-1|$ constraints |  Ref.
---|---|---
Solar-System bounds (Cassini mission) |  $\lesssim 0.2$ ∗ |  [62]
Supernova monopole radiation | $<10^{-2}$ |  [67]
Cluster density profiles (Max-BCG) | $<3.5\times 10^{-3}$ |  [68]
CMB spectrum | $<10^{-3}$ |  [69]
GW170817 (GW from BNS merger) | $<3\times 10^{-3}$ |  [43]
Period decay of binary system |  $\mathbf{<2.09\times 10^{-4}}$ † |  This work
|  |  (Tsujikawa model)
Cluster abundances | $<1.6\times 10^{-5}$ |  [70, 71]
CMB + BAO + $\sigma_{8}-\Omega_{m}$ relationship ∗∗ | $<3.7\times 10^{-6}$ |  [72]
Strong gravitational lensing (SLACS) | $<2.5\times 10^{-6}$ |  [73]
Redshift-space distortions | $<2.6\times 10^{-6}$ |  [74]
Distance indicators in dwarf galaxies | $<5\times 10^{-7}$ |  [75]
Table 3:  Comparison of the bounds on the cosmological scale amplitude
$f^{\prime}(R_{0})$ from different observations.
∗ This is obtained for the Tsujikawa model, when translated from the bound
$|f^{\prime}(R_{gal})-1|\lesssim 5\times 10^{-11}$ at the galactic scale.
† The Tsujikawa model gives the best constraint over other two models: the Hu-
Sawicki and the Starobinsky model.
∗∗ Taking into account cluster number counts (PSZ catalog) and weak-lensing
tomography measurements (CFHTLens). This analysis assumes the Hu-Sawicki
model.

Finally, the systems for observation that we considered here are quasi-stable
binary star systems. Observations of gravitational waves from binary merger
events by LIGO-Virgo provide the opportunity to test the gravity in the very
strong field regime and therefore can give more stringent bound on the $f(R)$
gravity. In our earlier work [43] we provided one such bound on $f(R)$ models
from the inspiral stage (dynamic but far away from the merger stage) of the
binary neutron star merger event GW170817. In the future observations, we hope
to see some detection of neutron star - black hole or neutron-star white-dwarf
merger events by LIGO- Virgo which will provide stronger bounds. Also it would
be interesting to see if our semi-classical approach can also be applied for
GW emission from inspiral or even in the merger stage.

## Acknowledgements

Research of SJ is partially supported by the SERB, DST, Govt. of India,
through a TARE fellowship grant no. TAR/2021/000354, hosted by the department
of Physics, Indian Institute of Technology Kharagpur.

## References

  * [1] S. Weinberg, _Gravitation and Cosmology, p.p:55-57_ , John Wiley and Sons, New York (1972).
  * [2] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Observation of gravitational waves from a binary black hole merger_ , _Phys. Rev. Lett._ 116 (2016) 061102.
  * [3] E.H.T. Collaboration, K. Akiyama, A. Alberdi, W. Alef, K. Asada, R. Azulay et al., _First m87 event horizon telescope results. i. the shadow of the supermassive black hole_ , _The Astrophysical Journal Letters_ 875 (2019) L1.
  * [4] A. De Felice and S. Tsujikawa, _f(R) theories_ , _Living Rev. Rel._ 13 (2010) 3 [1002.4928].
  * [5] S. Nojiri and S.D. Odintsov, _Unified cosmic history in modified gravity: From f(r) theory to lorentz non-invariant models_ , _Physics Reports_ 505 (2011) 59.
  * [6] S. Nojiri, S. Odintsov and V. Oikonomou, _Modified gravity theories on a nutshell: Inflation, bounce and late-time evolution_ , _Physics Reports_ 692 (2017) 1.
  * [7] S. Capozziello, _Curvature quintessence_ , _Int. J. Mod. Phys. D_ 11 (2002) 483 [gr-qc/0201033].
  * [8] S. Capozziello, V.F. Cardone, S. Carloni and A. Troisi, _Curvature quintessence matched with observational data_ , _Int. J. Mod. Phys. D_ 12 (2003) 1969 [astro-ph/0307018].
  * [9] S. Capozziello, S. Carloni and A. Troisi, _Quintessence without scalar fields_ , _Recent Res. Dev. Astron. Astrophys._ 1 (2003) 625 [astro-ph/0303041].
  * [10] S.M. Carroll, V. Duvvuri, M. Trodden and M.S. Turner, _Is cosmic speed-up due to new gravitational physics?_ , _Phys. Rev. D_ 70 (2004) 043528.
  * [11] A. Dolgov and M. Kawasaki, _Can modified gravity explain accelerated cosmic expansion?_ , _Physics Letters B_ 573 (2003) 1.
  * [12] V. Faraoni, _Matter instability in modified gravity_ , _Phys. Rev. D_ 74 (2006) 104017.
  * [13] I. Sawicki and W. Hu, _Stability of cosmological solutions in $f(r)$ models of gravity_, _Phys. Rev. D_ 75 (2007) 127502.
  * [14] Y.-S. Song, W. Hu and I. Sawicki, _Large scale structure of $f(r)$ gravity_, _Phys. Rev. D_ 75 (2007) 044004.
  * [15] T. Chiba, _1/r gravity and scalar-tensor gravity_ , _Physics Letters B_ 575 (2003) 1.
  * [16] G.J. Olmo, _The gravity lagrangian according to solar system experiments_ , _Phys. Rev. Lett._ 95 (2005) 261102.
  * [17] G.J. Olmo, _Post-newtonian constraints on $f(r)$ cosmologies in metric and palatini formalism_, _Phys. Rev. D_ 72 (2005) 083505.
  * [18] A.L. Erickcek, T.L. Smith and M. Kamionkowski, _Solar system tests do rule out $1/r$ gravity_, _Phys. Rev. D_ 74 (2006) 121501.
  * [19] S. Nojiri and S.D. Odintsov, _Modified gravity with negative and positive powers of curvature: Unification of inflation and cosmic acceleration_ , _Phys. Rev. D_ 68 (2003) 123512.
  * [20] S. Nojiri and S.D. Odintsov, _Modified $f(r)$ gravity consistent with realistic cosmology: From a matter dominated epoch to a dark energy universe_, _Phys. Rev. D_ 74 (2006) 086005.
  * [21] W. Hu and I. Sawicki, _Models of $f(r)$ cosmic acceleration that evade solar system tests_, _Phys. Rev. D_ 76 (2007) 064004.
  * [22] W. Hu and I. Sawicki, _Models of f(R) Cosmic Acceleration that Evade Solar-System Tests_ , _Phys. Rev._ D76 (2007) 064004 [0705.1158].
  * [23] A.A. Starobinsky, _Disappearing cosmological constant in f(R) gravity_ , _JETP Lett._ 86 (2007) 157 [0706.2041].
  * [24] S. Tsujikawa, _Observational signatures of $f(r)$ dark energy models that satisfy cosmological and local gravity constraints_, _Phys. Rev. D_ 77 (2008) 023507.
  * [25] S. Nojiri and S.D. Odintsov, _Modified $f(r)$ gravity unifying ${R}^{m}$ inflation with the $\Lambda\mathrm{CDM}$ epoch_, _Phys. Rev. D_ 77 (2008) 026007.
  * [26] G. Cognola, E. Elizalde, S. Nojiri, S.D. Odintsov, L. Sebastiani and S. Zerbini, _Class of viable modified $f(r)$ gravities describing inflation and the onset of accelerated expansion_, _Phys. Rev. D_ 77 (2008) 046009.
  * [27] S. Capozziello, C. Corda and M.F. De Laurentis, _Massive gravitational waves from f(r) theories of gravity: Potential detection with lisa_ , _Physics Letters B_ 669 (2008) 255.
  * [28] D. Liang, Y. Gong, S. Hou and Y. Liu, _Polarizations of gravitational waves in $f(r)$ gravity_, _Phys. Rev. D_ 95 (2017) 104034.
  * [29] C.P.L. Berry and J.R. Gair, _Linearized $f(r)$ gravity: Gravitational radiation and solar system tests_, _Phys. Rev. D_ 83 (2011) 104022.
  * [30] J. Khoury and A. Weltman, _Chameleon cosmology_ , _Phys. Rev. D_ 69 (2004) 044026.
  * [31] C. Burrage and J. Sakstein, _Tests of Chameleon Gravity_ , _Living Rev. Rel._ 21 (2018) 1 [1709.09071].
  * [32] X. Zhang, T. Liu and W. Zhao, _Gravitational radiation from compact binary systems in screened modified gravity_ , _Phys. Rev. D_ 95 (2017) 104027.
  * [33] T. Liu, X. Zhang, W. Zhao, K. Lin, C. Zhang, S. Zhang et al., _Waveforms of compact binary inspiral gravitational radiation in screened modified gravity_ , _Phys. Rev. D_ 98 (2018) 083023.
  * [34] M. Fierz, W.E. Pauli and P.A.M. Dirac, _On relativistic wave equations for particles of arbitrary spin in an electromagnetic field_ , _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 173 (1939) 211 [https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1939.0140].
  * [35] H. van Dam and M.J.G. Veltman, _Massive and massless Yang-Mills and gravitational fields_ , _Nucl. Phys. B_ 22 (1970) 397.
  * [36] V.I. Zakharov, _Linearized gravitation theory and the graviton mass_ , _JETP Lett._ 12 (1970) 312.
  * [37] K. Hinterbichler, _Theoretical Aspects of Massive Gravity_ , _Rev. Mod. Phys._ 84 (2012) 671 [1105.3735].
  * [38] C. de Rham, _Massive Gravity_ , _Living Rev. Rel._ 17 (2014) 7 [1401.4173].
  * [39] E. Mitsou, _Aspects of Infrared Non-local Modifications of General Relativity_ , Ph.D. thesis, Geneva U., 2015.  1504.04050.  10.1007/978-3-319-31729-8.
  * [40] A. Joyce, B. Jain, J. Khoury and M. Trodden, _Beyond the Cosmological Standard Model_ , _Phys. Rept._ 568 (2015) 1 [1407.0059].
  * [41] X. Zhang, T. Liu and W. Zhao, _Gravitational radiation from compact binary systems in screened modified gravity_ , _Phys. Rev. D_ 95 (2017) 104027.
  * [42] T. Liu, X. Zhang, W. Zhao, K. Lin, C. Zhang, S. Zhang et al., _Waveforms of compact binary inspiral gravitational radiation in screened modified gravity_ , _Phys. Rev. D_ 98 (2018) 083023.
  * [43] S. Jana and S. Mohanty, _Constraints on $f(r)$ theories of gravity from gw170817_, _Phys. Rev. D_ 99 (2019) 044056.
  * [44] R.A. Hulse and J.H. Taylor, _Discovery of a pulsar in a binary system_ , _Astrophys. J. Lett._ 195 (1975) L51.
  * [45] J.H. Taylor and J.M. Weisberg, _A new test of general relativity: Gravitational radiation and the binary pulsar PS R 1913+16_ , _Astrophys. J._ 253 (1982) 908.
  * [46] J.M. Weisberg and J.H. Taylor, _Observations of Post-Newtonian Timing Effects in the Binary Pulsar PSR 1913+16_ , _Phys. Rev. Lett._ 52 (1984) 1348.
  * [47] P.C. Peters and J. Mathews, _Gravitational radiation from point masses in a Keplerian orbit_ , _Phys. Rev._ 131 (1963) 435.
  * [48] J.M. Weisberg and Y. Huang, _Relativistic Measurements from Timing the Binary Pulsar PSR B1913+16_ , _Astrophys. J._ 829 (2016) 55 [1606.02744].
  * [49] M. Kramer et al., _Tests of general relativity from timing the double pulsar_ , _Science_ 314 (2006) 97 [astro-ph/0609417].
  * [50] J. Antoniadis et al., _A Massive Pulsar in a Compact Relativistic Binary_ , _Science_ 340 (2013) 6131 [1304.6875].
  * [51] P.C.C. Freire et al., _The relativistic pulsar-white dwarf binary PSR J1738+0333 II. The most stringent test of scalar-tensor gravity_ , _Mon. Not. Roy. Astron. Soc._ 423 (2012) 3328 [1205.1450].
  * [52] T. Damour and G. Esposito-Farèse, _Tensor-scalar gravity and binary-pulsar experiments_ , _Phys. Rev. D_ 54 (1996) 1474.
  * [53] T. Liu, X. Zhang and W. Zhao, _Constraining f(r) gravity in solar system, cosmology and binary pulsar systems_ , _Physics Letters B_ 777 (2018) 286.
  * [54] R.P. Feynman, _Feynman lectures on gravitation_ (1996).
  * [55] S. Weinberg, _Photons and Gravitons in $S$-Matrix Theory: Derivation of Charge Conservation and Equality of Gravitational and Inertial Mass_, _Phys. Rev._ 135 (1964) B1049.
  * [56] M.J.G. Veltman, _Quantum Theory of Gravitation_ , _Conf. Proc. C_ 7507281 (1975) 265.
  * [57] J.F. Donoghue, M.M. Ivanov and A. Shkerin, _EPFL Lectures on General Relativity as a Quantum Field Theory_ , 1702.00319.
  * [58] S. Mohanty and P. Kumar Panda, _Particle physics bounds from the Hulse-Taylor binary_ , _Phys. Rev. D_ 53 (1996) 5723 [hep-ph/9403205].
  * [59] S. Mohanty, _Astroparticle Physics and Cosmology: Perspectives in the Multimessenger Era_ , vol. 975 (2020), 10.1007/978-3-030-56201-4.
  * [60] T. Kumar Poddar, S. Mohanty and S. Jana, _Vector gauge boson radiation from compact binary systems in a gauged $L_{\mu}-L_{\tau}$ scenario_, _Phys. Rev. D_ 100 (2019) 123023 [1908.09732].
  * [61] T.K. Poddar, S. Mohanty and S. Jana, _Gravitational radiation from binary systems in massive graviton theories_ , _Journal of Cosmology and Astroparticle Physics_ 2022 (2022) 019.
  * [62] C.M. Will, _The Confrontation between General Relativity and Experiment_ , _Living Rev. Rel._ 17 (2014) 4 [1403.7377].
  * [63] C.M. Will, _The Confrontation between general relativity and experiment_ , _Living Rev. Rel._ 9 (2006) 3 [gr-qc/0510072].
  * [64] T.M. Tauris and T. Sennels, _Formation of the binary pulsars psr b2303+46 and psr j1141-6545 - young neutron stars with old white dwarf companions_ , _Astron. Astrophys._ 355 (2000) 236 [astro-ph/9909149].
  * [65] N.D.R. Bhat, M. Bailes and J.P.W. Verbiest, _Gravitational-radiation losses from the pulsar-white-dwarf binary PSR J1141-6545_ , _Phys. Rev. D_ 77 (2008) 124017 [0804.0956].
  * [66] J.P.W. Verbiest, N.D.R. Bhat and M. Bailes, _PSR J1141-6545: a powerful laboratory of GR and tensor-scalar theories of gravity_ , 10, 2012, DOI [1210.0224].
  * [67] A. Upadhye and J.H. Steffen, _Monopole radiation in modified gravity_ , 1306.6113.
  * [68] L. Lombriser, F. Schmidt, T. Baldauf, R. Mandelbaum, U.c.v. Seljak and R.E. Smith, _Cluster density profiles as a test of modified gravity_ , _Phys. Rev. D_ 85 (2012) 102001.
  * [69] Y.-S. Song, H. Peiris and W. Hu, _Cosmological constraints on $f(r)$ acceleration models_, _Phys. Rev. D_ 76 (2007) 063517.
  * [70] F. Schmidt, A. Vikhlinin and W. Hu, _Cluster constraints on $f(r)$ gravity_, _Phys. Rev. D_ 80 (2009) 083505.
  * [71] M. Cataneo, D. Rapetti, F. Schmidt, A.B. Mantz, S.W. Allen, D.E. Applegate et al., _New constraints on $f(r)$ gravity from clusters of galaxies_, _Phys. Rev. D_ 92 (2015) 044009.
  * [72] L. Boubekeur, E. Giusarma, O. Mena and H. Ramírez, _Current status of modified gravity_ , _Phys. Rev. D_ 90 (2014) 103512.
  * [73] T.L. Smith, _Testing gravity on kiloparsec scales with strong gravitational lenses_ , 0907.4829.
  * [74] L. Xu, _Constraints on $f(r)$ gravity through the redshift space distortion_, _Phys. Rev. D_ 91 (2015) 063008.
  * [75] B. Jain, V. Vikram and J. Sakstein, _Astrophysical tests of modified gravity: Constraints from distance indicators in the nearby universe_ , _The Astrophysical Journal_ 779 (2013) 39.