KCL-PH-TH/2020-70

CERN-TH-2020-198

# Constraints on extended gravity models through gravitational wave emission

Gaetano Lambiase  , Mairi Sakellariadou  , Antonio Stabile

(August 27, 2024)

###### Abstract

Using recent experimental results of detection of gravitational waves from the
binary black hole signals by Advanced LIGO and Advanced Virgo, we investigate
the propagation of gravitational waves in the context of fourth order gravity
nonminimally coupled to a massive scalar field. Gravitational radiation admits
extra massive modes of oscillation and we assume that the amplitude of these
modes is comparable to that of the massless mode. We derive the propagation
equation and effective mass for each degree of freedom and we infer, from the
current observational data, constraints on the free parameters of the gravity
models we considered. In particular, for $f(R)=R-R^{2}/R_{0}$, the constraint
obtained from the speed of gravitational waves is not compatible with the one
set by Solar System tests, which implies that amplitude of the massive modes
could not be detectable with current experiments on Earth.

##  1 Introduction

The recent observations of gravitational radiation [1, 2, 3, 4, 5, 6, 7, 8]
from binary black holes by the Advanced LIGO and the Advanced Virgo detectors,
offered the ultimate confirmation of Einstein’s theory of General Relativity
(GR) and a novel opportunity to test all models of extended gravity in a
strong field dynamical regime. These alternative models to GR have been widely
studied and tested to the regime in which the Newtonian and post-Newtonian
contributions describe very well the dynamics. In particular, they have been
constrained by Earth tests (e.g., gravitoelectromagnetism and torsion balance
experiments), Solar System tests and by astrophysical data on the variation of
the orbital period of binary systems. However, direct detection of
gravitational waves from coalescing binary black holes, in an extreme
gravitational field regime, offers a new test for extended theories of
gravity.

On August 2017, the association of gravitational wave GW170817 [9] and gamma-
ray burst GRB170817A [10] events allowed to verify experimentally, with
deviations smaller than a few $10^{-15}$ [11], that gravitational waves travel
at the speed of light $c$. This measurement has considerable relevance since
it provides new insight into fundamental physics and probes the speed of
gravitational waves over cosmological distances.

Einstein’s theory of GR has been tested to very high precision in the last
decade [12]. To date it is the best theory of gravitational interaction,
however despite its great success there are still open questions which make
the theory incomplete. Several observational data [13, 14, 15, 16, 17, 18]
probe that the Universe is spatially flat and is undergoing a period of
accelerated expansion. To understand and describe the current astrophysical
data two unrevealed ingredients are needed in order to achieve this
phenomenology: dark matter at galactic and extragalactic scales and dark
energy at cosmological scales. The extended gravity models have been
considered as a viable theoretical mechanism to explain cosmic acceleration
and galactic rotation curves. In such models, one extends only the geometric
sector, without introducing any exotic matter. In such models, gravity shows a
different behaviour either above (infrared modification) or below (ultraviolet
modification) a certain length scale, whilst the robust results of General
Relativity at local and Solar System scales are preserved [19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34] (see Refs. [35, 36] for a class of
models related to Horndeski theory).

Extended theories of gravity have been widely studied in the Newtonian, post-
Newtonian and Minkowskian limit [37, 38, 39, 40, 44, 46, 43, 41, 42, 45, 52,
47, 48, 49, 50, 51, 53, 54]. In general, one finds modifications of the
gravitational potential in Newtonian and post-Newtonian limit, whilst in
Minkowskian limit one obtains massive gravitational wave modes. Many studies,
in the weak-field limit, have been conducted in context on various
astrophysical systems as the galactic rotation curves, stellar hydrodynamics
and gravitational lensing. An interesting application concerns the study of
the Casimir effect [56, 55, 57, 58, 59, 60, 61].

Corrections to the gravitational Lagrangian were already considered by several
authors [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. From a
conceptual viewpoint, there is no reason _a priori_ to restrict the
gravitational Lagrangian to a linear function of the Ricci scalar minimally
coupled to matter. In particular, one may consider the generalization of
$f(R)$ models, where $R$ is the Ricci scalar, through generic functions
containing curvature invariants such as the _Ricci squared_
($R_{\alpha\beta}R^{\alpha\beta}$) or the _Riemann squared_
($R_{\alpha\beta\gamma\delta}R^{\alpha\beta\gamma\delta}$), which however are
not invariant due to the Gauss-Bonnet invariant
$R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\lambda\sigma}R^{\mu\nu\lambda\sigma}$.
Note that the same remark applies to the Weyl invariant
$C_{\alpha\beta\gamma\delta}C^{\alpha\beta\gamma\delta}$. Hence, one may add a
(massive propagating) scalar field coupled to geometry; this leads to the
scalar-tensor fourth order gravity.

Using recent experimental results of detection of gravitational waves from the
binary black hole signals by Advanced LIGO and Advanced Virgo, we will study
the propagation of gravitational waves in the context of fourth order gravity
nonminimally coupled to a massive scalar field, in the weak field limit (for
different models see Refs. [62, 63, 64, 65, 66]. In particular, we will impose
constraints on the free parameters of extended gravity models from the current
observational data.

The outline of this paper is the following. In Section 2, we give the action
of a scalar tensor fourth order gravity and write down the corresponding field
equations, which we then solve in the presence of matter within the weak-field
approximation. In Section 3, using the recent observations of gravitational
radiation from binary black holes by the Advanced LIGO and the Advanced Virgo
detectors we put the constraint on the free parameters of the teories. In
Section 4 ,we analyze some models of extended theory. We summarize our
conclusions in Section 5.

##  2 Scalar Tensor Fourth Order Gravity

Let us consider the action

$\displaystyle\mathcal{S}\,=\,\int
d^{4}x\sqrt{-g}\biggl{[}f(R,R_{\alpha\beta}R^{\alpha\beta},\phi)+\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}+\mathcal{X}\mathcal{L}_{\rm
m}\biggr{]}\leavevmode\nobreak\ ;$ (2.1)

$f$ is an unspecified function of the Ricci scalar $R$, the curvature
invariant $R_{\alpha\beta}R^{\alpha\beta}\,\equiv Y$ (where $R_{\mu\nu}$ is
the Ricci tensor), and the scalar field $\phi$. The $\mathcal{L}_{\rm m}$
denotes the minimally coupled ordinary matter Lagrangian density, $\omega$ is
a generic function of the scalar field $\phi$, and $g$ stands for the
determinant of the metric tensor $g_{\mu\nu}$ and $\mathcal{X}\,\equiv\,8\pi
G$. Note that unless otherwise specified, we use natural units such that $c=1$
and $\hslash=1$.

In the metric approach, the field equations are obtained by varying the action
(2.1) with respect to $g_{\mu\nu}$:

$\displaystyle
f_{R}R_{\mu\nu}-\frac{f+\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}}{2}g_{\mu\nu}-f_{R;\mu\nu}+g_{\mu\nu}\Box
f_{R}+2f_{Y}{R_{\mu}}^{\alpha}R_{\alpha\nu}$
$\displaystyle-2[f_{Y}{R^{\alpha}}_{(\mu}]_{;\nu)\alpha}+\Box[f_{Y}R_{\mu\nu}]+[f_{Y}R_{\alpha\beta}]^{;\alpha\beta}g_{\mu\nu}+\omega(\phi)\phi_{;\mu}\phi_{;\nu}$
$\displaystyle=$ $\displaystyle\mathcal{X}\,T_{\mu\nu}\leavevmode\nobreak\ ;$
(2.2)

$T_{\mu\nu}\,=\,-\frac{1}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{m})}{\delta
g^{\mu\nu}}$ is the the energy-momentum tensor of matter,
$f_{R}\,=\,\frac{\partial f}{\partial R}$, $f_{Y}\,=\,\frac{\partial
f}{\partial Y}$ and $\Box=\nabla_{\mu}\nabla^{\mu}={{}_{;\sigma}}^{;\sigma}$
is the D’Alembertian operator111The convention for the Ricci tensor is
$R_{\mu\nu}={R^{\sigma}}_{\mu\sigma\nu}$, while for the Riemann tensor we
define ${R^{\alpha}}_{\beta\mu\nu}=\Gamma^{\alpha}_{\beta\nu,\mu}+...$, with
the affinities being the usual Christoffel’s symbols of the metric:
$\Gamma^{\mu}_{\alpha\beta}=\frac{1}{2}g^{\mu\sigma}(g_{\alpha\sigma,\beta}+g_{\beta\sigma,\alpha}-g_{\alpha\beta,\sigma})$.
The adopted signature is $(+---)$.. The trace of the field equations (2)
reads:

$\displaystyle
f_{R}R+2f_{Y}R_{\alpha\beta}R^{\alpha\beta}-2f+\Box[3f_{R}+f_{Y}R]+2[f_{Y}R^{\alpha\beta}]_{;\alpha\beta}-\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}\,=\,\mathcal{X}\,T\leavevmode\nobreak\
;$ (2.3)

$T\,=\,T^{\sigma}_{\,\,\,\,\,\sigma}$ is the trace of energy-momentum tensor.
Varying the action (2.1) with respect to scalar field $\phi$ we get:

$\displaystyle
2\omega(\phi)\Box\phi+\omega_{\phi}(\phi)\phi_{;\alpha}\phi^{;\alpha}-f_{\phi}\,=\,0\leavevmode\nobreak\
,$ (2.4)

where $\omega_{\phi}(\phi)\,=\,\frac{d\omega(\phi)}{d\phi}$ and
$f_{\phi}\,=\,\frac{\partial f}{\partial\phi}$.

Let us analyse the field equations within the _weak-field approximation_ in a
Minkowski background $\eta_{\mu\nu}$:

$g_{\mu\nu}\,\sim\,\eta_{\mu\nu}+h_{\mu\nu}\leavevmode\nobreak\
\leavevmode\nobreak\ ,\leavevmode\nobreak\ \leavevmode\nobreak\
\phi\,\sim\,\phi^{(0)}+\varphi\leavevmode\nobreak\ .$

We develop the function $f$ as:

$\displaystyle f(R,R_{\alpha\beta}R^{\alpha\beta},\phi)\,\sim$ $\displaystyle
f_{R}(0,0,\phi^{(0)})\,R+\frac{f_{RR}(0,0,\phi^{(0)})}{2}\,R^{2}+\frac{f_{\phi\phi}(0,0,\phi^{(0)})}{2}\,(\phi-\phi^{(0)})^{2}$
(2.5)
$\displaystyle+f_{R\phi}(0,0,\phi^{(0)})R\,(\phi-\phi^{(0)})+f_{Y}(0,0,\phi^{(0)})R_{\alpha\beta}R^{\alpha\beta}\leavevmode\nobreak\
.$

Note that any other possible contribution to $f$ is negligible [37, 38, 39].
Using the harmonic gauge condition,
$g^{\rho\sigma}\Gamma^{\alpha}_{\,\,\,\rho\sigma}\,=\,0$, we get
$h_{\mu\sigma}^{\,\,\,\,\,\,\,,\sigma}-1/2\,h_{\,,\mu}\,=\,0$, and hence Eq.
(2) reads [45]:

$\displaystyle(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}h_{\mu\nu}-\biggl{[}\frac{{m_{R}}^{2}-{m_{Y}}^{2}}{3{m_{R}}^{2}}\,\partial^{2}_{\mu\nu}+\eta_{\mu\nu}\biggl{(}\frac{{m_{Y}}^{2}}{2}+\frac{{m_{R}}^{2}+2{m_{Y}}^{2}}{6{m_{R}}^{2}}\Box_{\eta}\biggr{)}\biggr{]}\Box_{\eta}h$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\,{m_{Y}}^{2}\,f_{R\phi}(0,0,\phi^{(0)})\,(\partial^{2}_{\mu\nu}-\eta_{\mu\nu}\Box_{\eta})\varphi\,=\,-2\,{m_{Y}}^{2}\,\mathcal{X}\,T_{\mu\nu}\leavevmode\nobreak\
,$ (2.6)
$\displaystyle(\Box_{\eta}+{m_{R}}^{2})\Box_{\eta}h+6{m_{R}}^{2}\,f_{R\phi}(0,0,\phi^{(0)})\,\Box_{\eta}\varphi\,=\,2{m_{R}}^{2}\,\mathcal{X}\,T\leavevmode\nobreak\
,$ (2.7)
$\displaystyle(\Box_{\eta}+{m_{\phi}}^{2})\varphi+\frac{f_{R\phi}(0,0,\phi^{(0)})}{2}\,\Box_{\eta}h\,=\,0\leavevmode\nobreak\
;$ (2.8)

$\Box_{\eta}$ is the D’Alambertian operator in flat space,
$h\,=\,{h^{\sigma}}_{\,\sigma}$ and we have set222One could define a new
gravitational constant:
$\mathcal{X}\,\rightarrow\,\mathcal{X}\,f_{R}(0,0,\phi^{(0)})$ and
$f_{R\phi}(0,0,\phi^{0})\,\rightarrow\,f_{R\phi}(0,0,\phi^{0})\,f_{R}(0,0,\phi^{(0)})$.
$f_{R}(0,0,\phi^{(0)})\,=\,1$ and $\omega(\phi^{(0)})\,=\,1/2$. The quantities
${m_{R}}^{2},{m_{Y}}^{2},{m_{\phi}}^{2}$ are defined333In our formalism we
consider the coefficients of $R^{2}$ and $\phi^{2}$ negative. For this reason,
inside of Eqs. (2.14) we have the minus signs. However, the definitions of the
masses Eqs. (2.14) are consistent with the formalism used by K. S. Stelle in
Ref. [67]. as follows:

$\displaystyle\begin{array}[]{ll}{m_{R}}^{2}\,\equiv-\frac{f_{R}(0,0,\phi^{(0)})}{3f_{RR}(0,0,\phi^{(0)})+2f_{Y}(0,0,\phi^{(0)})}\leavevmode\nobreak\
,\\\ \\\
{m_{Y}}^{2}\,\equiv\frac{f_{R}(0,0,\phi^{(0)})}{f_{Y}(0,0,\phi^{(0)})}\leavevmode\nobreak\
,\\\ \\\
{m_{\phi}}^{2}\,\equiv-\frac{f_{\phi\phi}(0,0,\phi^{(0)})}{2\omega(\phi^{(0)})}\leavevmode\nobreak\
.\end{array}$ (2.14)

The field equations (2), (2.7) and (2.8) generalise those of Ref. [40], since
in the latter there was no scalar field component. Let us also note that these
equations are the weak-field limit of the model discussed in Refs. [44, 45].

To solve Eqs. (2), (2.7) and (2.8) we introduce the auxiliary field
$\gamma_{\mu\nu}$, such that:

$\displaystyle(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}\gamma_{\mu\nu}=\,-2\,{m_{Y}}^{2}\,\mathcal{X}\,T_{\mu\nu}\leavevmode\nobreak\
.$

Then the fields $h_{\mu\nu}$ can be written as [45]:

$\displaystyle
h_{\mu\nu}=\gamma_{\mu\nu}-\frac{A(\partial)}{m^{2}_{Y}}\gamma-2\,{m_{R}}^{2}\,f_{R\phi}(0,0,\phi^{(0)})\,B(\partial)\,\varphi\leavevmode\nobreak\
,$ (2.15)

where

$\displaystyle A(\partial)$ $\displaystyle=$
$\displaystyle\biggl{[}\frac{{m_{R}}^{2}-{m_{Y}}^{2}}{3}\,\partial^{2}_{\mu\nu}+\eta_{\mu\nu}\biggl{(}\frac{{{m_{R}}^{2}\,m_{Y}}^{2}}{2}+\frac{{m_{R}}^{2}+2{m_{Y}}^{2}}{6}\Box_{\eta}\biggr{)}\biggr{]}(\Box_{\eta}+{m_{R}}^{2})^{-1}\,,$
$\displaystyle B(\partial)$ $\displaystyle=$
$\displaystyle\,\biggl{[}\Box_{\eta}^{-1}(\Box_{\eta}+{m_{R}}^{2})^{-1}\partial^{2}_{\mu\nu}+\frac{1}{2}\eta_{\mu\nu}(\Box_{\eta}+{m_{R}}^{2})^{-1}\biggl{]}.$

Using Eq. (2.15), Eqs. (2), (2.7) and (2.8) becomes:

$\displaystyle(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}\gamma_{\mu\nu}=\,-2\,{m_{Y}}^{2}\,\mathcal{X}\,T_{\mu\nu}\leavevmode\nobreak\
,$ (2.17)
$\displaystyle(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}\gamma=\,-2\,{m_{Y}}^{2}\,\mathcal{X}\,T\leavevmode\nobreak\
,$ (2.18)
$\displaystyle\frac{(\Box_{\eta}+{m_{\phi}}^{2})\varphi}{{m_{R}}^{2}f^{2}_{R\phi}(0,0,\phi^{(0)})}-3\,(\Box_{\eta}+{m_{R}}^{2})^{-1}\Box_{\eta}\varphi\,=\,-(\Box_{\eta}+{m_{R}}^{2})^{-1}\mathcal{X}\,T\,.$
(2.19)

Thus, Eqs. (2.7) and (2.8) have been decoupled. Let us now rewrite Eq. (2.19)
as:

$\displaystyle(\Box_{\eta}+m^{2}_{+})(\Box_{\eta}+m^{2}_{-})\varphi\,=\,-{m_{R}}^{2}f_{R\phi}(0,0,\phi^{(0)})\,\mathcal{X}\,T\leavevmode\nobreak\
;$ (2.20)
$m^{2}_{\pm}\equiv{m_{R}}^{2}\,w_{\pm}\bigl{(}\xi,\eta\bigl{)}\leavevmode\nobreak\
,$ (2.21)

where

$\displaystyle w_{\pm}(\xi,\eta)$ $\displaystyle=$
$\displaystyle\frac{1-\xi+\eta^{2}\pm\sqrt{(1-\xi+\eta^{2})^{2}-4\eta^{2}}}{2}\,,$
$\displaystyle\xi$ $\displaystyle=$ $\displaystyle
3{f_{R\phi}(0,0,\phi^{(0)})}^{2}\,,$ $\displaystyle\eta\,$ $\displaystyle=$
$\displaystyle\,\frac{m_{\phi}}{m_{R}}\,,$

with the constraint $\xi<1$.

Applying on the first member of Eq. (2.15) the operators:

$\displaystyle(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}\,\,\,\,\,\,\text{and}\,\,\,\,\,\,(\Box_{\eta}+m^{2}_{+})(\Box_{\eta}+m^{2}_{-})\,,$

using Eqs. (2.17), (2.18) and (2.20), we get

$\displaystyle(\Box_{\eta}+m^{2}_{+})(\Box_{\eta}+m^{2}_{-})(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}h_{\mu\nu}=-2\,{m_{Y}}^{2}\,\mathcal{X}(\Box_{\eta}+m^{2}_{+})(\Box_{\eta}+m^{2}_{-})\,T_{\mu\nu}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+2\,m^{2}_{R}\,\mathcal{X}\,A(\partial)(\Box_{\eta}+m^{2}_{+})(\Box_{\eta}+m^{2}_{-})T$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+2\,{m_{R}}^{4}\,f_{R\phi}(0,0,\phi^{(0)})^{2}\,\mathcal{X}\,B(\partial)(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}T\leavevmode\nobreak\
.$ (2.22)

Considering the associated homogeneous field equation:

$\displaystyle(\Box_{\eta}+m^{2}_{+})(\Box_{\eta}+m^{2}_{-})(\Box_{\eta}+{m_{Y}}^{2})\Box_{\eta}h_{\mu\nu}=0\,,$
(2.23)

we note that the dynamics of the perturbation $h_{\mu\nu}$ is made by four
different modes: a massles one, which is the standard graviton, and three
massive modes with masses $m_{+}$, $m_{-}$ and $m_{Y}$. These quantities are
free parameters which can be constrained via observational data. In our
analysis, we suppose that all masses are different.

##  3 Experimental constraints

For the graviton (massless mode) the dispersion relation is $E^{2}=p^{2}$,
while for the three massive modes, $E^{2}=p^{2}+m^{2}_{\alpha}$, with
$\alpha\equiv\\{+,-,Y\\}$, or equivalently:

$\displaystyle
E^{2}=p^{2}+\,m_{\alpha}^{2}\,,\,\,\,\,\,\,\alpha\equiv\\{+,-,Y\\}\,.$ (3.1)

The group velocity [68], $v=\frac{p}{E}$ associated with the dispersion
relation (3.1) is:

$\displaystyle
v^{g}_{\alpha}=\sqrt{1-\frac{m_{\alpha}^{2}}{E^{2}}}\,,\,\,\,\,\,\,\alpha\equiv\\{+,-,Y\\}\,,$
(3.2)

where $v$, $p$ and $E$ denote respectively, the velocity, momentum and energy.
Note that if the rest mass vanishes, $m_{\alpha}=0$, the group velocity is 1,
and we have a massless particle, i.e. the graviton. Whereas, for
$m_{\alpha}\neq 0$, the group velocity is less than that of the graviton,
since we have massive particle.

We will test the compatibility of Scalar Tensor Fourth Order Gravity theories
(2.1) using current experimental data. We assume that the amplitude of massive
modes is either comparable to that of the massless mode, or it lies below the
detection threshold. Thus, experiments should be able to detect such massive
modes independently of whether the massless mode has been detected. Note that
if the amplitude of the massive modes is smaller than that of the massless
mode, one would not expect a detection and consequently there would be no
constraints.

Using the recent experimental results on the detection of gravitational waves
(massless mode) by Advanced LIGO and Advanced Virgo [7, 69] and following the
results presented in the reference [69] we can impose the constraint on the
three massive modes :

$\displaystyle m_{\alpha}<\,0.583\times
10^{-22}\,\text{eV}\,,\,\,\,\,\,\,\alpha\equiv\\{+,-,Y\\}\,.$ (3.3)

For the next section, we indicate with $m_{0}$ the value of upper bound for
the graviton mass, i.e. $m_{0}\equiv 0.583\times
10^{-22}\,\text{eV}=2.96\times 10^{-16}\,\text{m}^{-1}$.

##  4 Scalar tensor fourth order gravity models

Case | EGM | Parameters
---|---|---
|  | $m^{2}_{R}$ | $m^{2}_{Y}$ | $m^{2}_{\phi}$ | $\xi$ | $\eta$ | $m^{2}_{\pm}$
A | $f(R)$ | $-\frac{f_{R}(0)}{3f_{RR}(0)}$ | $\infty$ | 0 | 0 | 0 | 0
B | $f(R,R_{\alpha\beta}R^{\alpha\beta})$ | $-\frac{f(0)}{3f_{RR}(0)+2f_{Y}(0)}$ | $\frac{f_{R}(0)}{f_{Y}(0)}$ | 0 | 0 | 0 | 0
C | $f(R,\phi)+\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}$ | $-\frac{f_{R}(0)}{3f_{RR}(0)}$ | $\infty$ | $-\frac{f_{\phi\phi}(0)}{2\omega(\phi^{(0)})}$ | $\frac{3{f_{R\phi}(0)}^{2}}{2\omega(\phi^{(0)})}$ | $\frac{m_{\phi}}{m_{R}}$ | $m^{2}_{R}w_{\pm}(\xi,\eta)$
D | $f(R,R_{\alpha\beta}R^{\alpha\beta},\phi)+\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}$ | $-\frac{f(0)}{3f_{RR}(0)+2f_{Y}(0)}$ | $\frac{f_{R}(0)}{f_{Y}(0)}$ | $-\frac{f_{\phi\phi}(0)}{2\omega(\phi^{(0)})}$ | $\frac{3{f_{R\phi}(0)}^{2}}{2\omega(\phi^{(0)})}$ | $\frac{m_{\phi}}{m_{R}}$ | $m^{2}_{R}w_{\pm}(\xi,\eta)$
Table 1:  Some models of extended theories of gravity.

Let us consider some extended gravity models studied in the literature:

  * •

Case A: $f(R)$ denotes a family of theories, each one defined by a different
function $f$ of the Ricci scalar $R$. The simplest case is for $f$ being equal
to Ricci scalar, obtaining just General Relativity. In this case the
characteristic scale (mass) $m_{R}$ (see, Table 1, case A) depends only on the
first and second derivatives of $f(R)$. Therefore, if we consider, for
example, a polynomial expression
$f(R)\,=\,R+\alpha\,R^{2}+\sum_{n\,=\,3}^{N}\alpha_{n}\,R^{n}$, only the terms
$R$ and $R^{2}$ make contributions. Using the constraint Eq. (3.3), i.e.
$m_{R}<m_{0}$, we get the following relation between the derivates of $f(R)$:

$-f_{RR}(0)>\frac{f_{R}(0)}{3\,m_{0}^{2}}.$ (4.1)

A particular case of an $f(R)$ theory is the Starobinsky model,
$f(R)\,=\,R-R^{2}/R_{0}$ [20], with $R_{0}$ a constant. Here the mass (2.14)
reads $m_{R}=\sqrt{R_{0}/6}$ ($m_{Y}=\infty,\,m_{\phi}=0$) and using Eq. (3.3)
we obtain:

$\displaystyle R_{0}<6\,m_{0}^{2}\approx\,5.3\times 10^{-31}\,\text{m}^{-2}.$
(4.2)

However, at present the best constraint on $R_{0}$ is obtained by Solar System
tests [70] and reads $R_{0}>2.5\times 10^{8}\,\text{m}^{-2}$; a constraint
with which (4.2) is not compatible. One may thus conclude that the amplitude
of massive modes is not detectable on Earth. Although the bounds we obtained
are weaker, they could still be of interest, if for instance there are
modifications of the effective form of $f(R)$ in different regions, through
the chameleon mechanism.

  * •

Case B: $f(R,\,R_{\alpha\beta}R^{\alpha\beta})$, namely we also include the
curvature invariant $R_{\alpha\beta}R^{\alpha\beta}$. In this case there are
two characteristic scales $m_{R}$ and $m_{Y}$. Using Eq. (3.3), $m_{R}<m_{0}$,
and $m_{Y}<m_{0}$, with the mass definition Eqs. (2.14), we obtain the
following relations between the derivatives:

$\displaystyle-
f_{RR}(0,0)>\frac{2}{3}f_{Y}(0,0)+\frac{f_{R}(0,0)}{3\,m_{0}^{2}}\ \ ,\ \
f_{Y}(0,0)>\frac{f_{R}(0,0)}{m_{0}^{2}}.$ (4.3)

As an illustration, let us consider
$f(R)\,=\,R-R^{2}/R_{0}+R_{\alpha\beta}R^{\alpha\beta}/R_{ic}$, where $R_{0}$
and $R_{ic}$ are constants. The masses (2.14) are then given by

$\displaystyle
m_{R}=\sqrt{\frac{R_{0}\,R_{ic}}{6R_{ic}-2R_{0}}}\,,\,\,\,\,m_{Y}=\sqrt{R_{ic}}\,\,\,\,\text{and}\,\,\,m_{\phi}=0\,.$

In this case the inequalities (4.3) become:

$\displaystyle R_{0}<\frac{6m_{0}^{2}\,R_{ic}}{2m_{0}^{2}+R_{ic}}\ \ ,\ \
R_{ic}<m_{0}^{2}\,,$

implying

$\displaystyle R_{0}<2\,m^{2}_{0}\approx\,1.8\times
10^{-31}\,\text{m}^{-2}\,,\,\,\,\,\,\,R_{ic}<m^{2}_{0}\approx 8.8\times
10^{-32}\,\text{m}^{-2}\,.$

This class of theories includes the case of Weyl squared, _i.e._
$C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}\,=\,2R_{\mu\nu}R^{\mu\nu}-\frac{2}{3}R^{2}$,
where there is only one characteristic scale since
$m_{R}\,\rightarrow\,\infty$.

  * •

Case C: Scalar-tensor models
$f(R,\,\phi)+\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}$. There there are two
effective scales $m_{+}$ and $m_{-}$ generated from the interaction between
gravity and the scalar field (see, Table 1, case C). Imposing Eq. (3.3),
$m_{+}<m_{0}$ and $m_{-}<m_{0}$, and after some algebra we get

$\displaystyle m_{R}\,m_{\phi}<m_{0}^{2}\ \ $ , $\displaystyle\ \
\bigg{(}1-\frac{m_{0}^{2}}{m_{R}^{2}}\bigg{)}\bigg{(}1-\frac{m_{\phi}^{2}}{m_{0}^{2}}\bigg{)}<\xi\leqslant\bigg{(}\frac{m_{\phi}}{m_{0}}-1\bigg{)}^{2}\,,$
(4.4) $\displaystyle m_{R}\,m_{\phi}<m_{0}^{2}\ \ $ , $\displaystyle\ \
0<\xi\leqslant\bigg{(}\frac{m_{\phi}}{m_{R}}-1\bigg{)}^{2}\,.$ (4.5)

Therefore, from Eqs. (2.14), (4.4) and (4.5) we find the following relations
for the derivatives of the unspecified function $f$:

$\displaystyle\frac{f_{RR}(0,\,\phi^{(0)})}{f_{\phi\phi}(0,\phi^{(0)})}>\,\frac{f_{R}(0,\,\phi^{(0)})}{6\,\omega(\phi^{(0)})\,m_{0}^{4}}\,,\,\,\,\,\,\,\,\,\Theta_{C}<3{f_{R\phi}^{2}(0,\,\phi^{(0)})}\leqslant\Lambda_{C}\,,$
(4.6)

where

$\displaystyle\Theta_{C}$ $\displaystyle\equiv$
$\displaystyle\bigg{(}1+\frac{3\,m_{0}^{2}\,f_{RR}(0,\,\phi^{(0)})}{f_{R}(0,\,\phi^{(0)})}\bigg{)}\bigg{(}1+\frac{f_{\phi\phi}(0,\phi^{(0)})}{2\omega(\phi^{(0)})\,m_{0}^{2}}\bigg{)}\,,$
$\displaystyle\Lambda_{C}$ $\displaystyle\equiv$
$\displaystyle\bigg{(}\sqrt{\frac{-f_{\phi\phi}(0,\phi^{(0)})}{2\,\omega(\phi^{(0)})\,m_{0}^{2}}}-1\bigg{)}^{2}\,.$

As an example, let us consider the the following scalar-tensor theory:

$\displaystyle f_{\rm
ST}(R,\phi)\,=R-\frac{R^{2}}{R_{0}}+\lambda\,\varphi\,R+\frac{1}{2}\varphi_{,\alpha}\varphi^{,\alpha}-\frac{1}{2}{\cal
M}^{2}\,\varphi^{2}\leavevmode\nobreak\ .$

In this case, the masses (2.14) are given by

$\displaystyle
m_{R}=\sqrt{\frac{R_{0}}{6}}\,,\,\,\,\,m_{Y}=\infty\,\,\,\,\text{and}\,\,\,m_{\phi}={\cal
M}\,.$

The relations Eqs. (4.6) read

$\displaystyle R_{0}<\frac{6\,m_{0}^{4}}{{\cal M}^{2}}\,,\ \
\bigg{(}1-\frac{6\,m^{2}_{0}}{R_{0}}\bigg{)}\bigg{(}1-\frac{{\cal
M}^{2}}{m_{0}^{2}}\bigg{)}<3\,\lambda^{2}\leqslant\bigg{(}\frac{{\cal
M}}{m_{0}}-1\bigg{)}^{2}\,.$

The first inequality shows a simple relation between the scalar curvature
$R_{0}$ and the mass of the scalar field ${\cal M}$.

  * •

Case D:
$f(R,\,R_{\alpha\beta}R^{\alpha\beta},\phi)+\omega(\phi)\phi_{;\alpha}\phi^{;\alpha}$,
for which we have three effective scales $m_{+}$ and $m_{-}$ and $m_{Y}$.
Imposing the constraints on all masses ($m_{+}<m_{0}$, $m_{-}<m_{0}$ and
$m_{Y}<m_{0}$), we get the following relations

$\displaystyle m_{R}\,m_{\phi}<m_{0}^{2}\,$ , $\displaystyle\ \
\bigg{(}1-\frac{m_{0}^{2}}{m_{R}^{2}}\bigg{)}\bigg{(}1-\frac{m_{\phi}^{2}}{m_{0}^{2}}\bigg{)}<\xi\leqslant\bigg{(}\frac{m_{\phi}}{m_{0}}-1\bigg{)}^{2}\,,$
$\displaystyle m_{R}\,m_{\phi}<m_{0}^{2}\,$ , $\displaystyle\ \
0<\xi\leqslant\bigg{(}\frac{m_{\phi}}{m_{R}}-1\bigg{)}^{2}\,,$ $\displaystyle\
\ m_{Y}<m_{0}$ .

Using the definitions of the masses (2.14), we get

$\displaystyle\frac{3\,f_{RR}(0,0,\,\phi^{(0)})+2\,f_{Y}(0,0,\,\phi^{(0)})}{f_{\phi\phi}(0,0,\phi^{(0)})}>\,\frac{f_{R}(0,0,\,\phi^{(0)})}{2\,\omega(\phi^{(0)})\,m_{0}^{4}}\,,$
$\displaystyle\Theta_{D}<3{f_{R\phi}^{2}(0,0,\,\phi^{(0)})}\leqslant\Lambda_{D}\,,$
(4.7) $\displaystyle
f_{Y}(0,0,\,\phi^{(0)})>\frac{f_{R}(0,0,\,\phi^{(0)})}{m_{0}^{2}}\,,$

where

$\displaystyle\Theta_{D}$ $\displaystyle\equiv$
$\displaystyle\bigg{(}1+\frac{m_{0}^{2}\,[3\,\,f_{RR}(0,0,\,\phi^{(0)})+2\,f_{Y}(0,0,\,\phi^{(0)})]}{f_{R}(0,0,\,\phi^{(0)})}\bigg{)}\bigg{(}1+\frac{f_{\phi\phi}(0,0,\phi^{(0)})}{2\omega(\phi^{(0)})\,m_{0}^{2}}\bigg{)}\,,$
$\displaystyle\Lambda_{D}$ $\displaystyle\equiv$
$\displaystyle\bigg{(}\sqrt{\frac{-f_{\phi\phi}(0,0,\phi^{(0)})}{2\,\omega(\phi^{(0)})}}-1\bigg{)}^{2}\,.$

In this case as example we consider the following model:

$\displaystyle
f(R,R_{\alpha\beta}R^{\alpha\beta},\phi)\,=R-\frac{R^{2}}{R_{0}}+\frac{R_{\alpha\beta}R^{\alpha\beta}}{R_{ic}}+\lambda\,\varphi\,R+\frac{1}{2}\varphi_{,\alpha}\varphi^{,\alpha}-\frac{1}{2}{\cal
M}^{2}\,\varphi^{2}\leavevmode\nobreak\ ,$

with masses

$\displaystyle
m_{R}=\sqrt{\frac{R_{0}\,R_{ic}}{6R_{ic}-2R_{0}}}\,,\,\,\,\,m_{Y}=\sqrt{R_{ic}}\,\,\,\,\text{and}\,\,\,m_{\phi}={\cal
M}\,.$

The relations (4.7) read

$\displaystyle
R_{0}\,<\frac{\,6\,m_{0}^{4}\,R_{ic}}{2\,m_{0}^{4}+R_{ic}\,{\cal M}^{2}}\,,$
$\displaystyle\bigg{(}1-2\frac{(3R_{ic}-R_{0})\,m^{2}_{0}}{R_{0}\,R_{ic}}\bigg{)}\bigg{(}1-\frac{{\cal
M}^{2}}{m_{0}^{2}}\bigg{)}<3\,\lambda^{2}\leqslant\bigg{(}\frac{{\cal
M}}{m_{0}}-1\bigg{)}^{2}\,,$ $\displaystyle R_{ic}<m_{0}^{2}\,,$

from which we derive the following inequality between the scalar curvature
$R_{0}$ and the mass of the scalar field ${\cal M}$:

$\displaystyle R_{0}\,<\frac{\,6\,m_{0}^{4}}{2\,m_{0}^{2}+{\cal M}^{2}}\,.$

In general, the constraints one can obtain are rather weak, see Eq. (4.2), but
can in principle be improved once further data of nearby binary black hole
signals are available.

##  5 Conclusions

We have studied the propagation of gravitational waves in the context of
fourth order gravity nonminimally coupled to a massive scalar field, in the
weak-field approximation.

We have seen that, in general, extended models of gravity have some massive
modes in addition to the massless mode (graviton) of General Relativity. In
our analysis, we assumed that these modes has the amplitude comparable to that
of the massless mode. Using the recent observations of gravitational radiation
from the binary black hole signals detected by Advanced LIGO and Advanced
Virgo, published in the catalog GWTC-1 [6], we have constrained the free
parameters of extended gravity models. These constraints are weaker compared
to the ones obtained from experimental data of a different origin [44].

For the $f(R)$-theories we have analysed, the constraint derived from
gravitational waves is not compatible with the one currently obtained from
Solar System experiments. Therefore, the amplitude of massive modes cannot be
detected on Earth. This could happen because, these massive modes are either
significantly suppressed by a small coupling parameter or by a large energy
scale that affects them.

Nevertheless, the analysis presented here is important for two reasons.
Firstly, there could exist a, yet unknown, screening mechanism operating on
Earth and Solar System scales, (for instance the chameleon screening [71, 72,
73, 74]), but could not manifest on larger (astrophysical) scales. In fact,
one may expect that the variation of speed of GWs might occur in high density
environments, but such a screening effect should be reduced over distances of
order 40 Mpc [63, 35]. Secondly, the approach we have followed here can be
used once further data of nearby binary black hole signals are available,
probably leading to stronger constraints. Moreover, further observations over
larger distances could provide limits on both screening mechanisms and higher
derivative corrections, in particular on the effective model used in the
analysis presented here.

## Acknowledgments

M.S. is supported in part by the Science and Technology Facility Council
(STFC), United Kingdom, under the research grant ST/P000258/1. The authors
would like to thank the anonymous referee for observations that improved the
manuscript.

## References

  * [1] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, 061102 (2016)
  * [2] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, 241103 (2016)
  * [3] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. X 6, 041015 (2016)
  * [4] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 118, 221101 (2017)
  * [5] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 119, 141101 (2017)
  * [6] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. X 9 (2019)
  * [7] GWTC-1 Documentation [LIGO Scientific and Virgo], https://doi.org/10.7935/82H3-HH23 (2018)
  * [8] R. Abbott et al. [LIGO Scientific and Virgo], [arXiv:2010.14527 [gr-qc]].
  * [9] B. P. Abbott et. al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 119 , 16 161101 (2017)
  * [10] A. Goldstein et. al., Astrophys. J. 848 (2017)
  * [11] B. P. Abbott et. al. [LIGO Scientific and Virgo], Astrophys. J. 848 (2017)
  * [12] C. M. Will, Living Rev. Rel. 17, 4 (2014)
  * [13] A. G. Riess et al., Astron. J. , 116, 1009-1038 ( 1998)
  * [14] S. Perlmutter et al., Astron. J. ,517, 565-586 (1999)
  * [15] S. Cole et al., Mon. Not. R. Astron. Soc, 362, 505-534 (2005)
  * [16] D. N. Spergel et al., ApJS, 170, 377-408 (2007)
  * [17] S. M. Carroll, W. H. Press, E. L. Turner, Annu. Rev. Astron. Astrophys., 30, 499-542 (1992)
  * [18] V. Sahni, A. Starobinski, Int. J. Mod. Phys, 9, 373-443 (2000)
  * [19] A.A. Starobinsky, Phys. Lett. B 91, 99 (1980)
  * [20] A. A. Starobinsky, Sov. Astron. Lett. 9, 302, (1983)
  * [21] Christian G. Boehmer, Tiberiu Harko, Francisco S.N. Lobo, Astropart. Phys., 29, 386 (2008)
  * [22] C.G. Boehmer, T. Harko, F.S.N. Lobo, JCAP 0803, 024 (2008)
  * [23] T.P. Sotiriou, V. Faraoni, Rev. Mod. Phys.  82, 451 (2010)
  * [24] A. De Felice, S. Tsujikawa, Living Rev. Rel. 13, 3 (2010)
  * [25] S. Nojiri, S.D. Odintsov, V.K. Oikonomou, Phys. Rept. 692, 1 (2017)
  * [26] A. Silvestri, M. Trodden, Rep. Prog. Phys. 72, 096901 (2009)
  * [27] R. Durrer, R. Maartens, Gen. Relat. Grav. 40, 301 (2008)
  * [28] M. Sami, Lect. Notes Phys. 720, 219 (2007)
  * [29] E.J. Copeland, M. Sami, Sh Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006)
  * [30] L. Amendola, S. Tsujikawa, Dark Energy: Theory and Observations (Cambridge University Press, Cambridge, 2010)
  * [31] G. Lambiase S. Mohanty, A Stabile Eur. Phys. J. C 78, 350 (2018).
  * [32] M. Zubair and F. Kousar, Eur. Phys. J. C 76 (2016) no.5, 254
  * [33] M. Zubair, F. Kousar and S. Bahamonde, Phys. Dark Univ. 14 (2016), 116-125
  * [34] M. Zubair, F. Kousar and S. Bahamonde, Eur. Phys. J. Plus 133 (2018) no.12, 523
  * [35] J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, Phys. Rev. Lett. 114, 211101 (2015)
  * [36] J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, J. Cosmol. Astropart. Phys. 02 (2015) 018
  * [37] A. Stabile, Phys. Rev. D. 82, 064021(12) (2010)
  * [38] A. Stabile, Phys. Rev. D, 82, 124026(7) (2010)
  * [39] A. Stabile, S. Capozziello, Phys. Rev. D , 87, 064002 (13) (2013)
  * [40] A. Stabile, S. Capozziello, arXiv: 1501.02187
  * [41] A. Stabile, An. Stabile, Phys. Rev. D , 85, 044014 (2012)
  * [42] A. Stabile, An. Stabile and S. Capozziello, Phys. Rev. D , 88, 124011(9) (2013)
  * [43] G. Lambiase, M. Sakellariadou, An. Stabile, JCAP 2013, 12, 020 (2013)
  * [44] S. Capozziello, G. Lambiase, M. Sakellariadou, A. Stabile, An. Stabile, Phys. Rev. D, 91, 044012 (15) (2015)
  * [45] G. Lambiase, M. Sakellariadou, A. Stabile, An. Stabile, JCAP 1507, 003 (2015)
  * [46] T. Biswas, E. Gerwick, T. Koivisto, A. Mazumdar, Phys. Rev. Lett. 108, 031101 (2012)
  * [47] L. Buoninfante, A. S. Koshelev, G. Lambiase and A. Mazumdar, JCAP 1809, 034 (2018)
  * [48] L. Buoninfante, A. S. Koshelev, G. Lambiase, J. Marto and A. Mazumdar, JCAP 1806, 014 (2018)
  * [49] L. Buoninfante, G. Lambiase, Y. Miyashita, W. Takebe and M. Yamaguchi, Phys. Rev. D 101, 084019 (2020)
  * [50] L. Buoninfante, A. S. Cornell, G. Harmsen, A. S. Koshelev, G. Lambiase, J. Marto and A. Mazumdar, Phys. Rev. D 98, 084041 (2018)
  * [51] L. Buoninfante, G. Lambiase and M. Yamaguchi, Phys. Rev. D 100, 026019 (2019)
  * [52] L. Buoninfante, G. Lambiase and A. Stabile, Eur. Phys. J. C 80, 122 (2020)
  * [53] G.M. Tino, L. Cacciapuoti, S. Capozziello, G. Lambiase, F. Sorrentino, Prog. Part. Nucl. Phys. 112, 103772 (2020)
  * [54] G. Lambiase , L. Mastrototaro, Astrophys. J. 904 (2020) 1, 19
  * [55] G. Lambiase, A. Stabile and An. Stabile, Phys. Rev. D 95, 084019 (2017)
  * [56] M. Blasone, G. Lambiase, L. Petruzziello and A. Stabile, Eur. Phys. J. C 78, no. 11, 976 (2018)
  * [57] L. Buoninfante, G. Lambiase, L. Petruzziello and A. Stabile, Eur. Phys. J. C 79, no. 1, 41 (2019)
  * [58] I.L. Buchbinder, P.M. Lavrov, S.D. Odintsov, Nucl.Phys. B 308 (1988) 191-202
  * [59] V.G. Bagrov, I.L. Buchbinder, S.D. Odintsov, Phys. Lett. B 184 (1987) 202-208
  * [60] S.D. Odintsov, Sov. Phys. J. 32 (1989) 458-461
  * [61] E. Elizalde, S.D. Odintsov, A. Romeo, A.A. Bytsenko, S. Zerbini, Singapore: World Scientific (1994) 319 p
  * [62] C. de Rham, A.J. Tolley, Phys. Rev. D 101, 063518 (2020)
  * [63] P. Creminelli, F. Vernizzi, Phys. Rev. Lett. 119, 251302 (2017)
  * [64] Y. Gong, E. Papantonopoulos, Z. Yi1, Eur. Phys. J. C 78, 738 (2018)
  * [65] T. Inagaki, M. Taniguchi, Int. J. Mod. Phys. D 29, 2050072 (2020)
  * [66] J. Sakstein, B. Jain, Phys. Rev. Lett. 119, 251303 (2017)
  * [67] K.S. Stelle, Phys.Rev.D 16 (1977) 953-969
  * [68] C. M. Will, Phys. Rev. D 57, 2061 (1998)
  * [69] B. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. D 100, 104036 (2019)
  * [70] C. P. L. Berry, J. R. Gair, Phys Rev. D  83, 104022 (2011)
  * [71] J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004); Phys. Rev. D 69, 044026 (2004)
  * [72] K. Hinterbichler and J. Khoury, Phys. Rev. Lett. 104, 231301 (2010)
  * [73] J. Sakstein, Phys. Rev. D 97, 064028 (2018). See also X. Zhang, W. Zhao, H. Huang, Y. Cai, Phys. Rev. D 93, 124003 (2016)
  * [74] P. Brax, C. van de Bruck, C. Davies, J. Khoury, A. Weltman, Phys. Rev D 70, 123518 (2004)