# A new $f(R)$ Gravity Model and properties of Gravitational Waves in it

Dhruba Jyoti Gogoi<EMAIL_ADDRESS>Department of Physics, Dibrugarh
University, Dibrugarh 786004, Assam, India  Umananda Dev Goswami
<EMAIL_ADDRESS>Department of Physics, Dibrugarh University, Dibrugarh
786004, Assam, India

###### Abstract

In this paper, we have introduced a new $f(R)$ gravity model as an attempt to
have a model with more parametric control, so that the model can be used to
explain the existing problems as well as to explore new directions in physics
of gravity, by properly constraining it with recent observational data. Here
basic aim is to study the properties of Gravitational Waves (GWs) in this new
model. In $f(R)$ gravity metric formalism, the model shows the existence of
scalar degree of freedom as like other $f(R)$ gravity models. Due to this
reason, there is a scalar mode of polarization of GWs present in the theory.
This polarization mode exists in a mixed state, of which one is transverse
massless breathing mode with non-vanishing trace and the other is massive
longitudinal mode. The longitudinal mode being massive, travels at speed less
than the usual tensor modes found in General Relativity (GR). Moreover, for a
better understanding of the model, we have studied the potential and mass of
scalar graviton in both Jordan frame and Einstein frame. This model can pass
the solar system tests and can explain primordial and present dark energy.
Also, we have put constraints on the model. It is found that the correlation
function for the third mode of polarization under certain mass scale predicted
by the model agrees well with the recent data of Pulsar Timing Arrays. It
seems that this new model would be useful in dealing with different existing
issues in the areas of astrophysics and cosmology.

Modified Gravity; Gravitational Waves; Constraints; Pulsar Timing Arrays

###### pacs:

04.30.Tv, 04.50.Kd

##  I Introduction

Einstein’s General Relativity (GR) has been the most widely accepted theory
capable of explaining a number of phenomena and the geometry of spacetime in
general. However, recent experimental observations showed the existence of
phenomena with deviations from GR predictions. Among them the present
acceleration of the universe has been a big problem in cosmology lacking of a
proper or satisfactory explanation till now. This problem was discovered
around 22 years ago with the help of type Ia supernovae observations
Riess_1998 ; Perlmutter1999 ; Bahcall1999 ; Kirshner_1999 . It put a big
question on the viability of GR, specially at cosmological scale. However, if
one still wishes to stick with GR then he has to bear with the problem of
mysterious invisible and exotic dark energy, which is responsible for around
76% of total energy content of the universe. Moreover, a theoretical problem
associated with GR is that it is not renormalizable based on the conventional
methods. In order to overcome these drawbacks of GR different modifications
have been proposed. In most of the cases the modifications were introduced to
solve some specific problems and as expected/ eventually they redefined the
spacetime geometry and imprinted some changes in other sectors/ranges also.
These changes can be a measure of how much the new theory is deviating from
GR. A very important result from GR is the Gravitational Waves (GWs). In
modified theories of gravity, the properties of GWs also change and can result
massive modes of polarization apart from the usual GR modes or tensor modes of
polarization corda01 . These massive modes travel with less speed than the
tensor modes. Besides the presence of massive modes of polarizations, the
generation and propagation of GWs are also affected in different modified
theories of gravity. The study of these variations can be a good tool to test
the modified theories of gravity. With the first detection of GWs in 2015 by
the LIGO and Virgo collaborations Abbott_2016  followed by many other
detections till now, a new and promising direction of studying gravitational
theories has begun. These experimental evidences of GWs put a new set of
constraints on GR as well as on modified theories of gravity.

A straight forward and simple modification to GR is the $f(R)$ theory of
gravity. In this theory, the Ricci scalar $R$ of Einstein - Hilbert action is
replaced by a function of $R$. Till now many models of $f(R)$ theory have been
proposed. Some successful models capable of explaining the drawbacks of GR
upto a significant range are the Starobinsky model Starobinsky_2007 , Hu-
Sawicki model Hu_2007 , and Tsujikawa model Tsujikawa_2008 ; Zhang_2006 ;
Cognola_2008 ; Linder_2009  or exponential gravity. The properties of GWs in
these models have been studied earlier. These studies show that GWs in $f(R)$
gravity vary significantly from those in GR. As mentioned above, in GR, there
are two polarization modes of GWs, viz., the tensor plus mode and tensor cross
mode. These modes are massless in nature and they propagate with the speed of
light in spacetime. In metric formalism $f(R)$ gravity, there exists a scalar
degree of freedom in the theory and hence the total degrees of freedom in the
theory increases, which leads to increase the polarization modes of GWs in
such theories corda01 ; corda02 . It is found that total number of
polarization modes of GWs that can exist in $f(R)$ theories of gravity is 3
Liang_2017 . Recent studies show that there exist a massless breathing mode
and a massive longitudinal mode in a mixed state as a single polarization mode
of GWs in $f(R)$ theories of gravity. The existence of this extra polarization
mode can be checked with the help of modified Newman-Penrose scalars and also
by geodesic deviation equations in the theory Liang_2017 . However, another
study shows that the existence of third mixed polarization mode of GWs are
model dependent and there exists a model in $f(R)$ gravity where the massive
longitudinal mode vanishes and the third polarization mode becomes a pure
massless breathing mode Gogoi_2019 .

In this work, we have used a new $f(R)$ gravity model as a toy model and
studied the behaviour of the potential and scalar field mass both in Jordan
frame and Einstein frame. We have also checked the polarization modes of GWs
present in this model and tried to constrain the model.

The paper is organized as follows. In the next section, we have introduced our
new $f(R)$ model along with motivations. The characteristics of this model and
the behaviour of the associated scalar field both in Jordan frame and Einstein
frame have been studied in this section. Also, in this section, we have
checked the viability of the model in terms of solar system tests. In the
third section, we have compared the model with two other most viable models
viz., the Starobinsky model and the Hu-Sawicki model. In the fourth section,
the model has been constrained. In section five, we have studied the
polarization modes of GWs present in the model by using the perturbed field
equation and the modified Newman-Penrose formalism. In the sixth section, we
have reviewed a way to detect the polarization modes of GWs experimentally and
discussed the possibilities of experimental validation of the model. In the
last section, we conclude the paper with a very brief discussion of the
results and the future aspects of the model in such type of studies.

##  II A New Model of f(R) Gravity

Although we have a pretty good number of $f(R)$ gravity models, no $f(R)$
gravity model can explain all cosmological and astrophysical aspects of the
present universe completely. Moreover, as like GR, different $f(R)$ gravity
models come with different types of drawbacks. However, we should mention that
including the Starobinsky and Hu-Sawicki models there are several viable
models in $f(R)$ gravity, which were proposed in order to overcome the
drawbacks of GR. Although the asymptotic behaviour of such viable models are
almost same, the functional forms are completely different. In such $f(R)$
gravity models the modifications of geometry is done in a unique manner by
different $f(R)$ functions, which might result different and unique
cosmological and astrophysical features in the intermediate curvature regime.
Thus a new $f(R)$ gravity model capable of explaining the drawbacks of GR
might have some different implications in different perspectives and also can
have drawbacks or anomalies in different realms. Furthermore, a $f(R)$ gravity
model with more parametric control is more suitable in this respect in the
sense that such a model can be constrained properly with the observational
data by controlling its parameters and hence can be used easily to overcome
the drawbacks of GR. With these motivations, here we introduce another $f(R)$
gravity model containing three parameters, over and above to the existing
ones, as given by:

$f(R)=R-\frac{\alpha}{\pi}\,R_{c}\cot^{-1}\\!\left(\tfrac{R_{c}^{2}}{R^{2}}\right)-\beta\,R_{c}\\!\left[1-\exp\left({-\,\tfrac{R}{R_{c}}}\right)\right],$
(1)

where $\alpha$ and $\beta$ are two dimensionless positive constants and
$R_{c}$ is a characteristic curvature constant having dimensions same as
curvature scalar $R$. This model has two correction terms:

$\frac{\alpha}{\pi}\,R_{c}\cot^{-1}\\!\left(\tfrac{R_{c}^{2}}{R^{2}}\right)\;\;\text{and}\;\;\beta\,R_{c}\\!\left[1-\exp\left({-\,\tfrac{R}{R_{c}}}\right)\right].$

The first correction factor is estimated by two parameters $\alpha$ and
$R_{c}$. Similarly, the second correction factor has also two parameters
$\beta$ and $R_{c}$ and it mimics the exponential $f(R)$ gravity model. We’ll
show that this model passes the basic requirements of a viable $f(R)$ gravity
model including the solar system tests.

The requirements for any $f(R)$ gravity model to describe the late-time dark
energy problem are Chen_2019 ; Amendola_2015 ; Amendola_2007 ; Bamba_2010 :

(1) A sufficient and suitable chameleon mechanism which allows $f(R)$ gravity
to pass the constraints of local systems. In case of our model, we’ll show
that it can pass the solar system tests. A detailed study of chameleon
mechanism in this model is beyond the scope of this paper.

(2) A late-time stable de-Sitter solution. The condition for the existence of
de-Sitter solution for a model is

$2f(R_{0})-R_{0}f^{\prime}(R_{0})=0,$ (2)

where $R_{0}$ is the de Sitter curvature. To ensure the stability of the de
Sitter solution, the model needs to satisfy the following condition motohashi
:

$\frac{f^{\prime}(R_{0})}{f^{\prime\prime}(R_{0})}>R_{0}.$ (3)

For simplicity, considering $R_{0}=R_{c}$ in the Eq. (2) and solving for
$\beta$ we find,

$\beta=-\frac{e(2\alpha-\pi(\alpha-2))}{(6-4e)\pi},$ (4)

where $e$ is the natural exponential. This is the de Sitter solution for the
case $R_{0}=R_{c}$. Now the stability condition (3) gives,

$(e\,\alpha+\pi\,\beta)(2\,e\,\alpha+2\,\pi\,\beta-e\,\pi)<0.$ (5)

Using the expression for $\beta$ from Eq. (4) in the above expression, we find
the allowed range of $\alpha$ as

$-1.68381<\alpha<0.367545.$ (6)

Thus for $R_{0}=R_{c}$, the model can have stable de Sitter solutions if
$\alpha$ lies in the above range. In section III, we explicitly showed that
the model has a stable de Sitter solution and also compared it with two other
viable models.

(3) A positive effective gravitational coupling, leading to $f^{\prime}(R)>0$.
By putting our model in this condition gives,

$1-\frac{2\,\alpha\,R\,R_{c}^{3}}{\pi R^{4}+\pi
R_{c}^{4}}-\beta\exp\left(-\,\tfrac{R}{R_{c}}\right)>0.$

(4) A stable cosmological perturbation and a positivity of the GWs for the
scalar mode, causing to $f^{\prime\prime}(R)>0$. Using our model in this
condition we find,

$\frac{\beta\exp\left(-\,\tfrac{R}{R_{c}}\right)}{R_{c}}-\frac{2\,\alpha\,R_{c}^{3}\left(R_{c}^{4}-3R^{4}\right)}{\pi\left(R^{4}+R_{c}^{4}\right)^{2}}>0.$

This condition ensures the absence of tachyonic instabilities in the model,
i.e. this condition results $m^{2}_{\phi}>0$.

(5) An asymptotic behaviour to the $\Lambda CDM$ model in the large curvature
region. In case of our model, at large curvature region i.e. at
$R\rightarrow\infty$, we have
$f(R)-R\rightarrow-\,\frac{1}{2}R_{c}(\alpha+2\beta)=\text{constant}$, which
mimicks the $\Lambda CDM$ model in the large curvature region. Again, at
$R\rightarrow 0$, we have $f(R)-R\rightarrow 0$.

Thus it is seen that this model is capable of explaining the late-time dark
energy problem. In the following part of this paper, we’ll study the behaviour
of scalar degrees of freedom of the model as well as the nature of GWs in it.

###  II.1 Scalar Degrees of Freedom in Jordan Frame

The action of a generic $f(R)$ gravity model is given as

$S=\dfrac{1}{2\kappa}\int d^{4}x\sqrt{-g}f(R)+\int
d^{4}x\sqrt{-g}\,\mathcal{L}_{m}\\!\left[g^{\mu\nu},\bar{\psi}\,\right]\\!.$
(7)

In the above equation, the function $f(R)$ stands for any arbitrary function
of Ricci curvature scalar R, $g_{\mu\nu}$ is the metric, $\kappa^{2}=8\pi
G=M_{pl}^{-2}$ and $\hbar=c=1$. Here $M_{pl}\approx 2\times
10^{18}\,\text{GeV}$ is the reduced Planck’s mass.
$\mathcal{L}_{m}\\!\left[g^{\mu\nu},\bar{\psi}\,\right]$ is the Lagrangian for
a matter field $\bar{\psi}$. Variation of the action (7) with respect to the
metric gives the following field equation:

$f^{\prime}(R)R_{\mu\nu}-\dfrac{1}{2}f(R)g_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}f^{\prime}(R)+g_{\mu\nu}\,\square
f^{\prime}(R)=\kappa^{2}\,T_{\mu\nu}(g^{\mu\nu},\bar{\psi}),$ (8)

where $\square\equiv\nabla^{\mu}\nabla_{\mu}$,
$T_{\mu\nu}(g^{\mu\nu},\bar{\psi})=\dfrac{-\,2}{\sqrt{-g}}\dfrac{\delta\left(\sqrt{-g}\,\mathcal{L}_{m}\\!\left[g^{\mu\nu},\bar{\psi}\,\right]\right)}{\delta
g^{\mu\nu}}$ is the matter energy-momentum tensor and
$f^{\prime}(R)=\partial_{R}f(R)$. Trace of Eq. (8) is

$f^{\prime}(R)R+3\,\square f^{\prime}(R)-2f(R)=\kappa^{2}\,T,$ (9)

where $T=g^{\mu\nu}T_{\mu\nu}$ is the trace of $T_{\mu\nu}$. It is seen that
the trace of the field equation is dynamical. This equation also indicates the
existence of an extra scalar degree of freedom in the theory. For a detailed
study about this degree of freedom we would like to use our model given in the
Eq. (1). Now, if we define a scalar field as

$\phi=f^{\prime}(R),$ (10)

then for our model the field $\phi$ becomes,

$\phi=1-\frac{2\alpha R_{c}^{3}}{\pi
R^{3}\\!\left(\frac{R_{c}^{4}}{R^{4}}+1\right)}-\beta\,\exp\left(-\,\tfrac{R}{R_{c}}\right).$
(11)

This shows that the scalar curvature $R$ can be expressed as a function of the
field $\phi$. From this definition of scalar field $\phi$, we may view the
trace Eq. (9) as an effective scalar field equation of Klein-Gordon with the
following identification Capo_2007 :

$\displaystyle\dfrac{dV}{d\phi}$
$\displaystyle\equiv\dfrac{1}{3}\Big{[}2f(R(\phi))-R(\phi)f^{\prime}(R(\phi))\Big{]}$
$\displaystyle\equiv\dfrac{1}{3}\Big{[}2f(R(\phi))-R(\phi)\phi\Big{]},$ (12)

where $V$ is the potential of the scalar field $\phi$. Thus, the trace Eq. (9)
can be written as a Klein-Gordon type equation for the scalar field $\phi$ as
given by,

$\square\phi=\dfrac{dV}{d\phi}+\dfrac{1}{3}\,\kappa^{2}\,T=\dfrac{dV_{eff}}{d\phi},$
(13)

where $V_{eff}$ is effective potential of the field and is define as

$\dfrac{dV_{eff}}{d\phi}=\dfrac{1}{3}\Big{[}2f(R(\phi))-R(\phi)\phi+\kappa^{2}\,T\Big{]}.$
(14)

At far away from the source or in absence of any matter source $V_{eff}\equiv
V$. Again, from the stationary condition:

$\dfrac{dV_{eff}}{d\phi}=0,$ (15)

we can have $\phi=\phi_{0}$ satisfying $\phi_{0}=f^{\prime}(R_{0})$. From this
condition, the mass of the scalar field (or the scalaron mass) can be obtained
by differentiating the Eq. (14) with respect to $\phi$ as

$m^{2}_{\phi}\equiv\dfrac{dV^{2}_{eff}}{d\phi^{2}}\Big{|}_{\phi=\phi_{0}}=\dfrac{1}{3}\\!\left[\dfrac{f^{\prime}(R_{0})}{f^{\prime\prime}(R_{0})}-R_{0}\right]\\!.$
(16)

Here $R_{0}$ is the background curvature corresponding to $\phi_{0}$. From the
above equation, we can see that avoiding of tachyonic instabilities demands
$\frac{f^{\prime}(R_{0})}{f^{\prime\prime}(R_{0})}-R_{0}\geq 0$ and to keep
the mass term finite we need $f^{\prime\prime}(R_{0})\neq 0$. For our model
the mass term is found as

$m^{2}_{\phi}=\left[\frac{R_{c}\exp\\!\left(R/R_{c}\right)\left(\pi\left(R^{4}+R_{c}^{4}\right)^{2}-8\alpha
R^{5}R_{c}^{3}\right)-\pi\beta(R+R_{c})\left(R^{4}+R_{c}^{4}\right)^{2}}{3\pi\beta\left(R^{4}+R_{c}^{4}\right)^{2}-6\alpha
R_{c}^{4}\exp\\!\left(R/R_{c}\right)\left(R_{c}^{4}-3R^{4}\right)}\right]_{R\,=\,R_{0}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!.$
(17)

For $R_{0}=0$, this equation gives

$m^{2}_{\phi}\big{|}_{R_{0}\,=\,0}=\frac{\pi(\beta-1)R_{c}}{6\alpha-3\pi\beta}.$
(18)

This shows that the mass of the scalar field is non-vanishing even at a large
distance away from the source or in the Minkowski space. The mass $m_{\phi}$
of the scalar field depends on the model parameters $\alpha$, $\beta$ and
$R_{c}$. Fig. 1 shows the variation of $m^{2}_{\phi}$ with respect to $R_{0}$
for different sets of model parameters. From the figure we see that, the mass
of the scalar field increases rapidly with the increasing value of the
background curvature after a hump in the curve for $0<R_{0}<1$ region, which
increases when the parameter $\alpha$ takes value closer to parameter $\beta$.
By increasing the difference between $\beta$ and $\alpha$ $($i.e. for
$\beta\gg\alpha)$ the hump can be minimized. An increase of $\alpha$ increases
the hump which occurs near the small curvature region as mentioned above and
comparatively decreases the mass of the scalar field at higher curvature
region.

Figure 1: Mass square of the scalar field as a function of the background
curvature for different sets of values of $\alpha$ and $\beta$ with
characteristic curvature constant $R_{c}=1$ in arbitrary units.

####  II.1.1 Scalar Tensor Equivalence of the Model

To see the origin of the scalar field in the theory, we would like to rewrite
the action (7) by introducing a new auxiliary scalar field $\psi$ as
Liang_2017 ; Teyssandier_1983

$S=\dfrac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\,\Big{[}f_{\psi}(\psi)R-\big{\\{}f_{\psi}(\psi)\psi-f(\psi)\big{\\}}\Big{]}+\int
d^{4}x\sqrt{-g}\,\mathcal{L}_{m}\\!\left[g^{\mu\nu},\bar{\psi}\;\right],$ (19)

where $f_{\psi}=\dfrac{\partial f\\!\left(\psi\right)}{\partial\psi}$. Now,
varying this equation with respect to the new auxiliary scalar field $\psi$ we
get,

$f_{\psi\psi}\\!\left(\psi\right)\left(R-\psi\right)=0.$ (20)

For finiteness of the previously defined mass square term of the scalar field
(see Eq. (16)), we have $f^{\prime\prime}\\!\left(R\right)\neq 0$, which is in
terms of $\psi$, $f_{\psi\psi}\\!\left(\psi\right)\neq 0$. With this condition
the above equation gives, $R=\psi$. Substituting of this result in the action
(19) we can recover the original action (7). Moreover, the quantum stability
condition demands that $f^{\prime\prime}\\!\left(R\right)\geq 0$. This along
with the finiteness condition of the mass of the scalar field demands that
$f^{\prime\prime}\\!\left(R\right)>0$. Thus, it is always possible to have a
scalar tensor representation of $f(R)$ theory of gravity. Redefining the
previously defined scalar field $\phi$ in terms of the new auxiliary field
$\psi$ as

$\phi=f_{\psi}\\!\left(\psi\right),$ (21)

the action (19) can be rewritten as

$S=\dfrac{1}{2\kappa^{2}}\int d^{4}x\sqrt{-g}\,\Big{[}\phi
R-U\left(\phi\right)\Big{]}+\int
d^{4}x\sqrt{-g}\,\mathcal{L}_{m}\\!\left[g_{\mu\nu},\bar{\psi}\;\right],$ (22)

where
$U\left(\phi\right)=f_{\psi}\\!\left(\psi\right)\psi-f\\!\left(\psi\right)=\phi\,\psi\\!\left(\phi\right)-f\\!\left(\psi\left(\phi\right)\right)$
is the potential of the scalar field. To be precise, this term
$f_{\psi}\\!\left(\psi\right)\psi-f\\!\left(\psi\right)$ originates the scalar
field. Unless and otherwise this term equals to zero, there exists a scalar
field in the theory. In terms of the Ricci scalar, this term reads,

$V=U\\!\left(\phi\right)=f^{\prime}\\!\left(R\right)R-f\\!\left(R\right)\big{|}_{R\,=\,R_{0}}.$
(23)

For $f\\!\left(R\right)=R$ and hence for $f^{\prime}\\!\left(R\right)=1$, the
action (22) recovers GR giving the potential term $V=0$. Using Eq. (1) in this
Eq. (23), the scalar field potential for our model can be obtained as

$V=\frac{\alpha}{\pi}\,R_{c}\cot^{-1}\\!\left(\tfrac{R_{c}^{2}}{R^{2}}\right)+R\left[-\,\frac{2\alpha
R_{c}^{3}}{\pi
R^{3}\left(\frac{R_{c}^{4}}{R^{4}}+1\right)}-\beta\exp\left(-\,\tfrac{R}{R_{c}}\right)+1\right]+\beta
R_{c}\left[1-\exp\left(-\tfrac{R}{R_{c}}\right)\right]-R\;\Big{|}_{R\,=\,R_{0}}.$
(24)

Figure 2: Variation of the potential (24) as a function of the background
curvature $R_{0}$ for different sets of values of $\alpha$ and $\beta$ with
$R_{c}=1$ in arbitrary units.

The variation of the potential $V$ with respect to $R_{0}$ for different
values of $\alpha$ and $\beta$ is shown in Fig. (2) in arbitrary units
considering $R_{c}=1$. The figure shows that the potential in Jordan frame
increases gradually with respect to the background curvature with some initial
deviations depending upon the values of $\alpha$ and $\beta$. In contrast to
the case of mass square of the scalar field, the potential shows some slight
amount of dip, but near to the same small curvature region, which increases
slowly when the value of $\alpha$ moves closer to the value of $\beta$. In
fact, this dip in the potential is responsible for the hump in the mass square
curve for the corresponding values of $\alpha$ and $\beta$. This dip region of
the potential curve almost eliminates in the case when $\alpha\ll\beta$.
Moreover, with the increasing values of both $\alpha$ and $\beta$, the
potential comparatively increases after the dip region or without the dip
region.

###  II.2 Model in Einstein Frame

In order see the behaviour of our model (1) in Einstein frame, which is
usually used to avoid the non - minimal coupling of gravity with the scalar
field, we would like to study the model in this frame also. In the Einstein
frame, the following conformal transformation of the metric is performed
Goswami2013 ; Chakraborty2019 :

$\tilde{g}_{\mu\nu}=f^{\prime}(R)g_{\mu\nu},$

which for our model takes the form:

$\tilde{g}_{\mu\nu}=\left[1-\frac{2\alpha R_{c}^{3}}{\pi
R^{3}\left(\frac{R_{c}^{4}}{R^{4}}+1\right)}-\beta\exp\left(-\,\tfrac{R}{R_{c}}\right)\right]\\!g_{\mu\nu}.$
(25)

Consequently, in this frame with $\mathcal{L}_{m}=0$, the action changes to
Goswami2013 ; Chakraborty2019

$S=\int
d^{4}x\sqrt{-g}\left[\dfrac{1}{2\kappa^{2}}\tilde{R}-\dfrac{1}{2}\,\tilde{g}^{\mu\nu}\,\nabla_{\mu}\phi_{E}\nabla_{\nu}\phi_{E}-V(\phi_{E})\right]\\!,$
(26)

where the scalar field

$\phi_{E}=-\,\sqrt{\dfrac{3}{2}}\dfrac{1}{\kappa}\ln
f^{\prime}(R)=-\,\sqrt{\dfrac{3}{2}}\dfrac{1}{\kappa}\ln\\!\left[1-\frac{2\alpha
R_{c}^{3}}{\pi
R^{3}\left(\frac{R_{c}^{4}}{R^{4}}+1\right)}-\beta\exp\left(-\,\tfrac{R}{R_{c}}\right)\right]$
(27)

and $V(\phi_{E})$ is potential of the scalar field in this frame, and is given
by,

$V(\phi_{E})=\dfrac{1}{2\kappa^{2}}\dfrac{U}{f^{\prime}(R)^{2}}=\dfrac{1}{2\kappa^{2}}\dfrac{f^{\prime}(R)R-f(R)}{f^{\prime}(R)^{2}}\Big{|}_{R=R_{0}}.$
(28)

The Eq. (27) shows the dependency of the Einstein frame scalar field
$\phi_{E}$ on the scalar curvature $R$. Using our model, this potential (28)
can be expressed as

$V(\phi_{E})=\dfrac{1}{2\kappa^{2}}\,\frac{\pi\,R_{c}\,\chi\,e^{x}\Big{[}-\pi\,\beta\,\chi\left(x+1\right)+e^{x}\\!\left(\pi\beta+\pi\beta\,x^{4}-2\,\alpha\,x^{2}\right)+\alpha\,\chi\,e^{x}\cot^{-1}\\!\left(x^{-\,2}\right)\Big{]}}{\Big{[}e^{x}\left(\pi
x^{4}-2\,\alpha\,x+\pi\right)-\pi\beta\chi\Big{]}^{2}},$ (29)

where $x=R_{0}/R_{c}$ and $\chi=x^{4}+1$. Hence, the mass square term of the
scalar field in the Einstein frame is

$\displaystyle m^{2}(\phi_{E})$
$\displaystyle=\dfrac{d^{2}V(\phi_{E})}{d\phi_{E}^{2}}=\dfrac{1}{3}\left[\dfrac{1}{f^{\prime\prime}(R)}+\dfrac{R}{f^{\prime}(R)}-\dfrac{4f(R)}{f^{\prime}(R)^{2}}\right]_{R\,=\,R_{0}}$
$\displaystyle=\frac{1}{3}R_{c}\left[\frac{\pi\chi^{2}e^{x}}{2\,\alpha\,e^{x}\left(3x^{4}-1\right)+\pi\beta\chi^{2}}+\frac{\pi\chi
e^{x}x}{e^{x}\left(\pi x^{4}-2\,\alpha\,x+\pi\right)-\pi\beta\chi}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{4R_{c}\left[-\,\pi^{-1}\alpha\cot^{-1}\left(x^{-2}\right)+\beta\left(e^{-x}-1\right)+x\right]}{3\left[2\,\alpha\,x\left(\pi
x^{4}+\pi\right)^{-1}+\beta e^{-x}-1\right]^{2}}.$ (30)

In Minkowski space i.e. for $R_{0}=0$, this equation takes the form:

$m^{2}(\phi_{E})=\frac{\pi R_{c}}{3(\pi\beta-2\,\alpha)}.$ (31)

Figure 3: Variation of scalar field potential as a function of background
curvature in Einstein frame for different values of $\alpha$ and $\beta$
parameters with $R_{c}=1$ in arbitrary units.

Figure 4: Variation of scalaron mass square as a function of background
curvature in Einstein frame for different values of $\alpha$ and $\beta$
parameters with $R_{c}=1$ in arbitrary units.

Figure 5: Variation of scalar field $\phi_{E}$ as a function of scalar
curvature in Einstein frame for different values of $\alpha$ and $\beta$
parameters with $R_{c}=1$ in arbitrary units.

Figure 6: Variation of scalar field potential as a function of field
$\phi_{E}$ in Einstein frame for different values of $\alpha$ and $\beta$
parameters with $R_{c}=1$ in arbitrary units.

From Eq.s (29) and (II.2), we see that although the expressions for the scalar
field potential and scalaron mass square in Einstein frame are little bit
complicated in comparison to that in the Jordan frame, their variations as a
function of $R_{0}$ are found to be almost similar as seen from Fig. 3 and
Fig. 4 respectively. That is, the behaviours of potential and the mass term of
the scalar field are almost identical in both frames. In Fig. 5, the variation
of field $\phi_{E}$ as a function of background curvature is shown. The field
$\phi_{E}$ at $R_{0}=0$ and $R_{0}\rightarrow\infty$ is independent of
$\alpha$, whereas it is independent of $\beta$ only at
$R_{0}\rightarrow\infty$. Moreover, $\phi_{E}$ is non-zero at $R_{0}=0$ and
tends to zero at $R_{0}\rightarrow\infty$, which is obvious from it’s
expression. Again from Eq.s (11) and (27) it is clear that all these
behaviours of $\phi_{E}$ should be applicable to $\phi$ also, but with
positive values of $\phi$ for all values of $R_{0}$. Variation of the
potential (29) as a function of $\phi_{E}$ is shown in Fig. 6. It is seen from
this figure that the minimum of the potential moves towards the higher value
of $\phi_{E}$ when $\alpha\ll\beta$ than the case when $\alpha\sim\beta$.
Obviously, similar behaviour can be attributed to the potential (24) as a
function of the field $\phi$. Thus, because of the similarity of behaviours of
the scalar field, it’s potential and mass square term in both Jordan and
Einstein frames, the rest of the study in this paper is done in the Jordan
frame only.

###  II.3 Solar System Tests of the Model

It is possible to recover GR by introducing the Chameleon mechanism in the
theory. In this mechanism, the scalar field $\phi=f^{\prime}(R)$ is coupled
with the matter density of the environment. Thus, when a model is used inside
the solar system, due to presence of matter density, the scalar field coupled
with it gains mass and hence allows the model to pass the solar system tests.
Clearly, this mechanism implies that the functional form $f(R)$ should have a
very closer value to the Ricci scalar $R$, for $R$ above or equal to the solar
system scale. A model is considered viable and consistent if it passes the
solar system tests. Guo has introduced several methods to test whether an
$f(R)$ gravity model passes the solar system tests or not in Jordan frame
Guo_2014 . According to Guo, a model can pass solar system tests if it
satisfies the following conditions:

$\left|\,\dfrac{f(R)-R}{R}\,\right|\ll 1,$ (32)
$\left|\,f^{\prime}(R)-1\,\right|\ll 1,$ (33) $Rf^{\prime\prime}(R)\ll 1.$
(34)

We’ve calculated the above functions numerically for our model (see Table 1).
These functions are plotted against background curvature in the units of
$R_{c}$ for different parameters (see Fig. 7). These indicate that the model
can be made to pass the solar system tests by increasing the ratio
$R_{0}/R_{c}$ or by simply decreasing the parameter $R_{c}$. However, a simple
and effective way to make the model solar system viable is to decrease all the
parameters (i.e. $\alpha,\beta$ and $R_{c}$) sufficiently (see Table 1). Thus,
within a viable range of parameters, the model can easily pass the solar
system tests. This is another advantage of this model, which allows us to
enlist the model as a solar system viable model.

Figure 7: Plots of solar system test functions with respect to $R_{0}/R_{c}$. Table 1: Numerical values of solar system test functions for specific values of $R_{0}/R_{c},\alpha$ and $\beta$. Values of Model parameters |  $\left|\dfrac{f(R)-R}{R}\right|$ |  $\left|f^{\prime}(R)-1\right|$ |  $Rf^{\prime\prime}(R)$
---|---|---|---
$R_{0}/R_{c}=0.00148,\alpha=0.1500,\beta=0.500$ | 0.49970 | 0.49940 | 0.00060
$R_{0}/R_{c}=0.05000,\alpha=0.0450,\beta=0.500$ | 0.48842 | 0.47705 | 0.02235
$R_{0}/R_{c}=0.50000,\alpha=0.1500,\beta=0.500$ | 0.41686 | 0.34820 | 0.11727
$R_{0}/R_{c}=1.00000,\alpha=0.0500,\beta=0.080$ | 0.06307 | 0.04534 | 0.04534
$R_{0}/R_{c}=1.50000,\alpha=0.0050,\beta=0.044$ | 0.02401 | 0.01061 | 0.01657
$R_{0}/R_{c}=2.00000,\alpha=0.0005,\beta=0.008$ | 0.00356 | 0.00112 | 0.00226

##  III Comparison of the Model with other viable Models

In this section, we would like to compare the model with two other viable
models viz., the Starobinsky model and the Hu- Sawicki model in terms of
stability and nature of the models in local regime. At first, we would like
study the de Sitter stability of our toy model. For our model the Eq. (2)
takes the form:

$\frac{2\,\alpha\,x^{2}}{x^{4}+1}+\pi\left[\beta
e^{-x}(x+2)+x-2\beta\right]=2\alpha\cot^{-1}\\!\left(x^{-\,2}\right)\\!,$

where $x=R_{0}/R_{c}$ and $R_{c}\neq 0$. On solving this equation for $\beta$,
we get,

$\beta=-\,\frac{e^{x}\left(-\,\pi x^{5}-2\alpha
x^{2}+2\alpha\cot^{-1}\left(\frac{1}{x^{2}}\right)+2\alpha
x^{4}\cot^{-1}\left(\frac{1}{x^{2}}\right)-\pi
x\right)}{\pi\left(-\,x+2e^{x}-2\right)\left(x^{4}+1\right)}.$ (35)

The contour plot of $\beta$ as a function of $x$ and $\alpha$ is shown in Fig.
8.

Figure 8: $\beta$ as a function of $x$ and $\alpha$. The contour shows the de
Sitter solutions of the model.

The Starobinsky model and the Hu-Sawicki model are defined respectively as

$f_{S}(R)=R-s_{1}m_{s}^{2}\left[1-\left(\frac{R^{2}}{m_{s}^{4}}+1\right)^{-\,n}\right],$
(36)

and

$f_{H}(R)=R-\frac{m^{2}c_{1}\left(\frac{R}{m^{2}}\right)^{\mu}}{c_{2}\left(\frac{R}{m^{2}}\right)^{\mu}+1}.$
(37)

The de Sitter solution for the Starobinsky model is

$s_{1}=\frac{x_{s}\left(x_{s}^{2}+1\right)^{n+1}}{2\left(x_{s}^{2}\left(x_{s}^{2}+1\right)^{n}+\left(x_{s}^{2}+1\right)^{n}-nx_{s}^{2}-x_{s}^{2}-1\right)},$
(38)

where $x_{s}=R_{0}/m_{s}^{2}$. Similarly, for the Hu-Sawicki model, using
$\mu=1$ the de Sitter solution is found to be,

$c_{2}=\frac{c_{1}+\sqrt{(c_{1}-1)c_{1}}-1}{x_{h}}.$ (39)

The contour plots of the parameters $s_{1}$ and $c_{2}$ respectively for the
Starobinsky model and the Hu-Sawicki model are shown in Fig. 9. From the
contours it can be seen that for very small values of $n$ and $x_{s}$ in the
Starobinsky model, we do not have de Sitter solutions. The parameter $s_{1}$
rises very rapidly for the increase of the parameters $n$ and $x_{s}$. For the
Hu-Sawicki model, we see that, there is no de Sitter solutions present for
$c_{1}<1$ and the parameter $c_{2}$ rises very slowly for increments of
$c_{1}$ beyond $1$. On the other hand, for our toy model, we see that $\alpha$
and $x$ have a wider region of de Sitter solutions. The other parameter
$\beta$ has a higher value than $\alpha$ in the de Sitter solution space
towards higher values of $x$.

Figure 9: Figure on left side shows the de Sitter solutions of the Starobinsky
model and figure on right side shows the de Sitter solutions for the Hu-
Sawicki model.

The de Sitter space defined by the Eq. (35) is stable if it satisfies the
condition (3). Our model (1) in Eq. (3) leads to the inequality:

$\left[2\alpha
e^{x}\left(3x^{4}-1\right)+\pi\beta\left(x^{4}+1\right)^{2}\right]\left[\pi\beta(x+1)\left(x^{4}+1\right)^{2}+e^{x}\\{8\alpha
x^{5}-\pi\left(x^{4}+1\right)^{2}\\}\right]<0.$ (40)

By using Eq. (35) in the above inequality, we can further reduce it to the
following form:

$\displaystyle\left(x-2e^{x}+2\right)^{2}\left[\pi
x\left(x^{4}+1\right)^{2}-2\alpha\left(x^{4}+1\right)^{2}\cot^{-1}\left(\frac{1}{x^{2}}\right)+2\alpha
A\right]$
$\displaystyle\left[2\alpha(x+1)\left(x^{4}+1\right)^{2}\cot^{-1}\left(\frac{1}{x^{2}}\right)+\pi\left(-x^{2}-2x+2e^{x}-2\right)\left(x^{4}+1\right)^{2}-2\alpha
Bx^{2}\right]>0,$ (41)

where
$A=\left(x^{6}-3x^{5}-6x^{4}+e^{x}\left(6x^{4}-2\right)+x^{2}+x+2\right)$ and
$B=\left(x^{5}-3x^{4}+8\left(e^{x}-1\right)x^{3}+x+1\right)$. When the model
satisfies this Eq. (III), it can have stable de Sitter solutions. The
stability region of the model in parameter space is shown in Fig.10. Again, in
order to get oscillatory behaviour around the de Sitter space, the following
condition needs to be satisfied motohashi :

$\frac{f^{\prime}(R_{0})}{f^{\prime\prime}(R_{0})}>\dfrac{25R_{0}}{16}.$ (42)

Using our model in this inequality, we get,

$\frac{\left(x^{4}+1\right)\\{e^{x}\left(\pi x^{4}-2\alpha
x+\pi\right)-\pi\beta\left(x^{4}+1\right)\\}}{2\alpha
e^{x}\left(3x^{4}-1\right)+\pi\beta\left(x^{4}+1\right)^{2}}>\frac{25x}{16}$
(43)

Now, eliminating $\beta$ by using Eq. (35), we may rewrite the above
expression as

$-\,\frac{2\left(x^{4}+1\right)\\{-\pi\left(-x+e^{x}-1\right)\left(x^{4}+1\right)-\alpha\left(x^{4}+1\right)\cot^{-1}\left(\frac{1}{x^{2}}\right)+2\alpha\left(e^{x}-1\right)x\\}}{\pi
x\left(x^{4}+1\right)^{2}-2\alpha\left(x^{4}+1\right)^{2}\cot^{-1}\left(\frac{1}{x^{2}}\right)+2\alpha\\{x^{6}-3x^{5}-6x^{4}+e^{x}\left(6x^{4}-2\right)+x^{2}+x+2\\}}>\frac{25x}{16}.$
(44)

If the toy model satisfies the above condition, it will have a stable de
Sitter solution as well as oscillatory behaviour around the de Sitter
solution. Similarly, for the Starobinsky model, the stability condition (3)
takes the form:

$s_{1}n\\{(2n+1)x_{s}^{2}-1\\}\\{4s_{1}n(n+1)x_{s}^{3}-\left(x_{s}^{2}+1\right)^{n+2}\\}<0.$
(45)

Using Eq. (38) in the above inequality, we obtain,

$C\left[x^{4}\\{-2n^{2}+\left(x_{s}^{2}+1\right)^{n}-3n-1\\}+\left(x_{s}^{2}+1\right)^{n}+x_{s}^{2}\\{2\left(x_{s}^{2}+1\right)^{n}-n-2\\}-1\right]>0,$
(46)

where
$C=nx_{s}\left(x_{s}^{2}+1\right)^{2n+2}\\{(2n+1)x_{s}^{2}-1\\}\left[\left(x_{s}^{2}+1\right)^{n}+x_{s}^{2}\\{\left(x_{s}^{2}+1\right)^{n}-n-1\\}-1\right]^{2}$.
The Starobinsky model gives stable de Sitter solutions for any parameter set
satisfying the above inequality. For this model, the condition (42) gives,

$C\left[x_{s}^{4}\\{-50n^{2}+16\left(x_{s}^{2}+1\right)^{n}-16-57n\\}+x_{s}^{2}\\{32\left(x_{s}^{2}+1\right)^{n}-32-7n\\}+16\left(x_{s}^{2}+1\right)^{n}-16\right]>0.$
(47)

Any parameter sets in the Starobinsky model satisfying this condition will
give oscillations around de Sitter solutions in de Sitter space. Finally for
the Hu-Sawicki model, condition (3) gives:

$c_{1}\mu
x_{h}^{\mu+1}\\{c_{2}(\mu+1)x_{h}^{\mu}-\mu+1\\}\left[x_{h}\left(c_{2}x_{h}^{\mu}+1\right){}^{3}-c_{1}\mu
x_{h}^{\mu}\\{c_{2}(\mu+2)x_{h}^{\mu}-\mu+2\\}\right]>0.$ (48)

Now using Eq. (39) in the above expression and considering $\mu=1$, we get,

$\left(c_{1}-1\right)x_{h}\\{\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}+\left(c_{1}-1\right)x_{h}\\}\\{2\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}+2c_{1}x_{h}-x_{h}\\}>0.$
(49)

Similarly, for the Hu-Sawicki model, the condition for existence of
oscillations around stable de Sitter solutions gives:

$c_{1}^{2}x_{h}\\{\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}+\left(c_{1}-1\right)x_{h}\\}\\{\left(32c_{1}-41\right)\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}+\left(32c_{1}^{2}-57c_{1}+25\right)x_{h}\\}>0.$
(50)

We have plotted these inequalities to check the stability region of the models
in parameter space in in Figs. 11 and 12.

Figure 10: Stability region of the model (1). Figure on the left side shows
the stability region allowed by Eq. (III) and figure on the right side shows
the stability region allowed by Eq. (44).

Figure 11: Stability region of the Starobinsky model. Figure on the left side
shows the stability region allowed by Eq. (46) and figure on the right side
shows the stability region allowed by Eq. (47).

Figure 12: Stability region of the Hu-Sawicki model. Figure on the left side
shows the stability region allowed by Eq. (49) and figure on the right side
shows the stability region allowed by Eq. (50).

From the Fig. 10, it is seen that the toy model has a small region of
instability for lower values of $x$ and comparatively higher values of
$\alpha$ (greater than 0.4 approximately). For higher values of $x$, the model
has a continuous stability region in the parameter space. In case of the
Starobinsky model, from Fig. 11 it is seen that there is a small region of
instability in the parameter space. But in case of Hu-Sawicki model, the model
is stable only for the values of $c_{1}>1$ (Fig. 12). This suggests that, the
toy model has a wider range of stability region in the parameter space. Since,
the model also has a comparatively larger region admitting Eq. (42), it can be
easily concluded that the toy model with a selected set of parameters can pass
the stability conditions as well as the condition for having oscillatory
solutions in the de Sitter space.

Now, we would like to compare the bahviour of the toy model with the
Starobinsky and Hu-Sawicki models using the functions used in the expressions
(32), (33) and (34). These expressions ensure the viability of the models in
local systems. Although, in the previous section, we have chosen the
parameters freely to check the solar system viabilities, now we shall study
the same and compare with the other two models in de Sitter sitter stability
regions of the parameter spaces. However, it should be noted that most of the
parameter sets used in that sections are found to belong automatically in the
de Sitter sitter stability region of the model. The function used in relation
(32) takes the following respective form in the toy model, Starobinsky model
and Hu-Sawicki model in the de Sitter solution regime,

$\dfrac{f(R_{0})-R_{0}}{R_{0}}=\frac{\alpha\left(x^{4}+1\right)\cot^{-1}\left(\frac{1}{x^{2}}\right)-\left(e^{x}-1\right)\left(\pi
x^{4}+2\alpha x+\pi\right)}{\pi\left(-x+2e^{x}-2\right)\left(x^{4}+1\right)},$
(51)

$\dfrac{f_{S}(R_{0})-R_{0}}{R_{0}}=-\frac{\left(x_{s}^{2}+1\right)\left(\left(x_{s}^{2}+1\right){}^{n}-1\right)}{2\left(\left(x_{s}^{2}+1\right){}^{n}+x_{s}^{2}\left(\left(x_{s}^{2}+1\right){}^{n}-n-1\right)-1\right)}$
(52)

and

$\dfrac{f_{H}(R_{0})-R_{0}}{R_{0}}=-\frac{c_{1}\left(x_{h}\right)_{h}^{\mu+1}}{\sqrt{\left(c_{1}-1\right)c_{1}\left(x_{h}\right)_{h}^{2}}\left(x_{h}\right)_{h}^{\mu}+\left(c_{1}-1\right)\left(x_{h}\right)_{h}^{\mu+1}+\left(x_{h}\right)_{h}^{2}}.$
(53)

The second function in (33) gives,

$f^{\prime}(R_{0})-1=\frac{2\alpha\left(x^{4}+1\right)\cot^{-1}\left(\frac{1}{x^{2}}\right)-x\left(\pi\left(x^{4}+1\right)+4\alpha\left(e^{x}-1\right)\right)}{\pi\left(-x+2e^{x}-2\right)\left(x^{4}+1\right)},$
(54)

$f_{S}^{\prime}(R_{0})-1=\frac{nx_{s}^{2}}{-\left(x_{s}^{2}+1\right){}^{n}+x_{s}^{2}\left(-\left(x_{s}^{2}+1\right){}^{n}+n+1\right)+1}$
(55)

and

$f_{H}^{\prime}(R_{0})-1=-\frac{c_{1}\mu
x_{h}^{\mu+3}}{\left(\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}x_{h}^{\mu}+\left(c_{1}-1\right)x_{h}^{\mu+1}+x_{h}^{2}\right){}^{2}}.$
(56)

The third function for the respective models are calculated as:

$R_{0}f^{\prime\prime}(R_{0})=\frac{2\alpha
x\left(3x^{4}-1\right)}{\pi\left(x^{4}+1\right)^{2}}-\frac{x\left(-\pi
x^{5}-2\alpha x^{2}+2\alpha\cot^{-1}\left(\frac{1}{x^{2}}\right)+2\alpha
x^{4}\cot^{-1}\left(\frac{1}{x^{2}}\right)-\pi
x\right)}{\pi\left(-x+2e^{x}-2\right)\left(x^{4}+1\right)},$ (57)

$R_{0}f_{S}^{\prime\prime}(R_{0})=\frac{nx_{s}^{2}\left((2n+1)x_{s}^{2}-1\right)}{\left(x_{s}^{2}+1\right)\left(\left(x_{s}^{2}+1\right){}^{n}+x_{s}^{2}\left(\left(x_{s}^{2}+1\right){}^{n}-n-1\right)-1\right)}$
(58)

and

$R_{0}f_{H}^{\prime\prime}(R_{0})=\frac{c_{1}\mu
x_{h}^{\mu+3}\left((\mu+1)\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}x_{h}^{\mu}+\left(c_{1}-1\right)(\mu+1)x_{h}^{\mu+1}-(\mu-1)x_{h}^{2}\right)}{\left(\sqrt{\left(c_{1}-1\right)c_{1}x_{h}^{2}}x_{h}^{\mu}+\left(c_{1}-1\right)x_{h}^{\mu+1}+x_{h}^{2}\right){}^{3}}.$
(59)

These functions for the said models are compared in Fig. 13. It is seen that
the models are capable of passing the solar system tests in the de Sitter
stable regime for higher values of $x$. The Hu-Sawicki model in this regime
shows a constant behaviour. However, the Starobinsky model and the toy model
give the higher values of the test functions for lower $x$ values. In this
case, the behaviour of the toy model still indicates the ability to pass the
solar system tests. One interesting point to note that for the toy model and
the Starobinsky model, the asymptotic behaviours of the test functions are
almost same. These results suggest that the behaviour of the toy model is
closer to the Starobinsky model. However, to make a clear conclusion, a
detailed study of the models is required and this is beyond the scope of this
manuscript.

Figure 13: Solar system test functions with respect to $x$ for the toy model
(with $\alpha=0.3$), Starobinsky model (with $n$ = 2 at the upper panel and
$n=1$ at the lower panel) and Hu-Sawicki model (with $\mu=1$ and $c_{1}=1.5$
at the upper panel and $c_{1}=1.1$ at the lower panel).

##  IV Constraints on the Model

For the viability of a model in $f(R)$ gravity it is utmost necessary to
impose constraints on the model on the basis of different available
observational data. A model which passes such constraints are considered as a
viable model in $f(R)$ gravity. Starobinsky model and Hu Sawicki model are
examples of two viable models in this context as mentioned earlier. There are
several ways to constrain an $f(R)$ gravity model Boubekeur_2014 ; Desai_2018
; Jana_2019 ; Gupta_2019 ; Chiba_2007 ; Cataneo_2015 ; Nunes_2017 . A
constrained model is helpful to study different implications of the model. In
this section, we will try to constrain our toy model using the results
published in Boubekeur_2014 ; Desai_2018  and Jana_2019 . In Ref.
Boubekeur_2014 , authors carried out a Markov chain Monte Carlo (MCMC)
analysis for GWs from Hu Sawicki model using the data sets of cosmic microwave
background (CMB) and baryon acoustic oscillations (BAO) together with the
independent constraints on the relationship between the matter clustering
amplitude $\sigma_{8}$ and the matter mass-energy density $\Omega_{m}$ from
Planck Sunyaev-Zeldovich (PSZ) cluster number counts and also from the
CFHTLens weak lensing tomography measurements. Combining CMB, BAO and
$\sigma_{8}-\Omega_{m}$ relationship from the PSZ catalog Ade_2014 , they
obtained a bound which is still better than the bounds obtained from the GW
event GW170817 Jana_2019 . The bound on the parameter $f^{\prime}(R)$ reported
by them is

$-3.7\times 10^{-6}<f^{\prime}(R)-1<3.7\times 10^{-6},$ (60)

with $95\%$ confidence level at upper bound Boubekeur_2014 . On the other hand
in Ref. Desai_2018  a constraint was introduced on the Compton wavelength
$\lambda_{g}$ of the graviton. From their study, we have a constraint on
$\lambda_{g}^{-1}$ as given by,

$0~{}\text{m}^{-1}<\lambda_{g}^{-1}<1.098901099\times
10^{-23}~{}\text{m}^{-1}$ (61)

with $90\%$ confidence level on upper bound Desai_2018 .

Now, we have computed the values of $f^{\prime}(R)-1$ and $\lambda_{g}^{-1}$
for our model by taking into consideration of above cited respective upper
bounds with the corresponding confidence levels to constraint our model
parameters $R_{c},\alpha$ and $\beta$. The results of this computation along
with the contour plots are shown in Fig. 14. Here we have not considered the
constraint on the model parameters coming from the fact that at higher
curvatures it goes near to the cosmological constant. However, the parameters
are chosen from the stability region of the parameter space. It is seen that
the model can be a viable one within a proper range of variables. The figure
shows the contours with $95\%$ confidence level for $f^{\prime}(R)-1$ and
$90\%$ confidence level for $\lambda_{g}^{-1}$ (the larger contour) and with
$68\%$ confidence level for the both (the smaller contour). The central point
denotes the boundary value for both the parameters $f^{\prime}(R)-1$ and
$\lambda_{g}^{-1}$ and any value lower than the boundary value is viable. This
point corresponds to the galaxy cluster Abell 1689 data Desai_2018 . In the
plots we have considered three sets of the parameters and we see that the
smaller values of $\alpha$ allow the model to pass the constraints easily. In
the Fig. 14, we have also shown 3 other points corresponding to galaxy
clusters Abell 262, Abell 1991 and Abell 383 data from the Ref. Gupta_2019 .
All these points lie within the confidence level contours.

Figure 14: Contours with $95\%$ confidence level on the upper bounds of
$f^{\prime}(R)-1$ Boubekeur_2014  and $90\%$ confidence level on the upper
bounds of $\lambda_{g}^{-1}$ Desai_2018  (the larger contour) and with $68\%$
confidence level on the both (the smaller contour). The central red dot
denotes the $\lambda_{g}$ corresponding to galaxy cluster Abell 1689
Desai_2018 , blue dot corresponds to Abell 262, black one corresponds to Abell
1991 and yellow one corresponds to Abell 383 Gupta_2019 . data

The model can be constrained by using the GWs event GW170817 also. In a recent
study Jana_2019 , $f(R)$ gravity was constrained by using the GW170817. They
provided a bound on $f^{\prime}(R)$, which is

$-3\times 10^{-3}<f^{\prime}(R)-1<3\times 10^{-3}.$ (62)

Using our toy model in this expression, we find,

$-\,\frac{3}{1000}<\frac{2\alpha x}{\pi x^{4}+\pi}+\beta
e^{-x}<\frac{3}{1000}.$ (63)

Now, using Eq. (35) in the above expression, we can have,

$-\,\frac{3}{1000}<\frac{\pi
x\left(x^{4}+1\right)-2\alpha\left(x^{4}+1\right)\cot^{-1}\left(\frac{1}{x^{2}}\right)+4\alpha\left(e^{x}-1\right)x}{\pi\left(-x+2e^{x}-2\right)\left(x^{4}+1\right)}<\frac{3}{1000}.$
(64)

Choosing $x=7.5$ (an arbitrary point in the stability region from Fig. 10) and
considering $\alpha$ to be a positive quantity, the above equation reduces to,

$0\leq\alpha<0.743783.$ (65)

Figure 15: Plot on left shows the variation of $\lambda_{g}$ with respect to
$R_{0}$, plot in middle shows the variation of $f^{\prime}(R)-1$ with respect
to $R_{0}$ and plot on right shows variation of $m_{g}$ with respect to
$f^{\prime}(R)-1$ for (i) $\alpha=0.150,\beta=0.500,R_{c}=30\times
10^{-50}~{}m^{-2}$; (ii) $\alpha=0.045,\beta=0.500,R_{c}=70\times
10^{-50}~{}m^{-2}$; (iii) $\alpha=0.300,\beta=0.850,R_{c}=10\times
10^{-50}~{}m^{-2}.$ These are the constrained set of parameters of our model
that is done on the basis of upper bounds on $f^{\prime}(R)-1$ Boubekeur_2014
and $\lambda_{g}$ Desai_2018

.

In Fig. 15, we have shown the variations of $\lambda_{g}$ with respect to
$R_{0}$, $f^{\prime}(R)-1$ with respect to $R_{0}$ and $m_{g}$ with respect to
$f^{\prime}(R)-1$ for the values of the model parameters used in the contour
plots in Fig. 14. The hump in the mass of the scalar field encountered in the
Fig. 1 and in the Fig. 6 are also present in the $\lambda_{g}$ vs. $R_{0}$
curves. However, as mentioned earlier, this hump vanishes when
$(\beta-\alpha)\gg 0$. For higher curvatures, $\lambda_{g}$ rapidly moves
towards zero. The function $f^{\prime}(R)-1$ also has significantly higher
values near the Minkowski spacetime and as soon as the background Ricci
curvature increases, this model dependent function decreases rapidly. This
nature of the model is suitable for overcoming the local system constraints.
Again as seen from the figure, with the increase of the function
$f^{\prime}(R)-1$ mass of the scalar field decreases initially at a faster
rate and then becomes almost constant at later stage.

##  V Polarization Modes of GWs in the Model

In this section, we wish to check the polarization modes of GWs in the model.
In presence of massive polarization mode, it would be easy to constraint the
model using the experimental results. To explore the polarization modes of GWs
in the model, at first we’ll introduce the perturbation to the field equation.

###  V.1 Perturbation to the Field Equation

If there are propagating GWs in spacetime, then they perturbs the metric
around its background value. Considering the background metric as
$\bar{g}_{\mu\nu}$ we may express the spacetime metric to the first order of
perturbation value $h_{\mu\nu}$, which is usually usually very small, as

$g_{\mu\nu}=\bar{g}_{\mu\nu}+h_{\mu\nu},\;\;\mbox{where}\;\left|\,h_{\mu\nu}\,\right|\ll\left|\,\bar{g}_{\mu\nu}\,\right|.$
(66)

Now, expanding the Ricci tensor and the Ricci scalar upto the first order of
$h_{\mu\nu}$, we may write:

$\displaystyle R_{\mu\nu}$ $\displaystyle\simeq\tilde{R}_{\mu\nu}+\delta
R_{\mu\nu}+\mathcal{O}(h^{2})$
$\displaystyle=\tilde{R}_{\mu\nu}-\dfrac{1}{2}(\nabla_{\mu}\nabla_{\nu}h-\nabla_{\mu}\nabla^{\lambda}h_{\lambda\nu}-\nabla_{\nu}\nabla^{\lambda}h_{\mu\lambda}+\square
h_{\mu\nu})+\mathcal{O}(h^{2})$ (67)

and

$\displaystyle R$ $\displaystyle\simeq\tilde{R}+\delta R+\mathcal{O}(h^{2})$
$\displaystyle=\tilde{R}-\square
h+\nabla^{\mu}\nabla^{\nu}h_{\mu\nu}-\tilde{R}_{\mu\nu}h^{\mu\nu}+\mathcal{O}(h^{2}).$
(68)

Similarly, we may write for the $f(R)$ and $f^{\prime}(R)$ as

$\displaystyle f(R)$ $\displaystyle\simeq
f(\tilde{R})+f^{\prime}(\tilde{R})\,\delta R+\mathcal{O}(h^{4}),$ (69)
$\displaystyle f^{\prime}(R)$ $\displaystyle\simeq
f^{\prime}(\tilde{R})+f^{\prime\prime}(\tilde{R})\,\delta
R+\mathcal{O}(h^{4}),$ (70)

where $\tilde{R}$ is some constant curvature. Thus, due to the perturbation in
spacetime the trace equation (9) can be rewritten as

$3f^{\prime\prime}(\tilde{R})\,\square\delta
R+\left[f^{\prime\prime}(\tilde{R})\tilde{R}-f^{\prime}(\tilde{R})\right]\\!\delta
R=0,$ (71)

where we have used $T_{\mu\nu}=0$ for the empty space or far away from the
source. Fixing the gauge to be harmonic gauge with
$\nabla_{\mu}h^{\mu}_{\nu}=\dfrac{1}{2}\,\nabla_{\nu}h$, which after operating
by $\nabla^{\nu}$ we find,

$\nabla^{\mu}\nabla^{\nu}h_{\mu\nu}=\dfrac{1}{2}\,\square h.$ (72)

An important point to be mentioned here is that, the Eq. (8) is also satisfied
by another solution: $R_{\mu\nu}=\Lambda g_{\mu\nu}=\tilde{R}_{\mu\nu}$,
giving

$2f^{\prime}(\tilde{R})\tilde{R}_{\mu\nu}-g_{\mu\nu}f(\tilde{R})=\kappa^{2}\,T_{\mu\nu}.$
(73)

This equation actually corresponds to the Eq. (15), which is the stationary
condition used earlier in the Sec. II. In empty space, this equation has the
form:

$2f^{\prime}(\tilde{R})\tilde{R}_{\mu\nu}-g_{\mu\nu}f(\tilde{R})=0$ (74)

and it leads to have the equation,

$2f(\tilde{R})-\tilde{R}f^{\prime}(\tilde{R})=0.$ (75)

This is the stationary condition (15) of the scalar field potential in the
empty space corresponding to the constant curvature $\tilde{R}$ of spacetime.
Using the Eq.s (68), (72) and (75) in Eq. (71), we get

$3f^{\prime\prime}(\tilde{R})\,\square^{2}h+\left(\frac{5f(\tilde{R})f^{\prime\prime}(\tilde{R})}{f^{\prime}(\tilde{R})}-f^{\prime}(\tilde{R})\right)\\!\square
h+\left(\frac{2f(\tilde{R})^{2}f^{\prime\prime}(\tilde{R})}{f^{\prime}(\tilde{R})^{2}}-f(\tilde{R})\right)\\!h=0.$
(76)

Now, we would like to define $\square h=m^{2}h$, where $m$ is the mass of the
associated scalar field. Using this definition in the above equation we
obtain,

$3f^{\prime\prime}(\tilde{R})\,m^{4}+\left(\frac{5f(\tilde{R})f^{\prime\prime}(\tilde{R})}{f^{\prime}(\tilde{R})}-f^{\prime}(\tilde{R})\right)\\!m^{2}+\left(\frac{2f(\tilde{R})^{2}f^{\prime\prime}(\tilde{R})}{f^{\prime}(\tilde{R})^{2}}-f(\tilde{R})\right)=0.$
(77)

This is a quadratic equation in $m^{2}$ and solution for $m^{2}$ gives,

$m^{2}=\frac{f^{\prime}(\tilde{R})}{3f^{\prime\prime}(\tilde{R})}-\frac{2f(\tilde{R})}{3f^{\prime}(\tilde{R})}=\dfrac{1}{3}\left[\dfrac{f^{\prime}(\tilde{R})}{f^{\prime\prime}(\tilde{R})}-\tilde{R}\right]\\!,$
(78)

and

$m^{2}=-\,\frac{f(\tilde{R})}{f^{\prime}(\tilde{R})}=-\,\dfrac{\tilde{R}}{2}.$
(79)

We see that the second solution corresponds to tachyonic scalar field which
becomes zero in the Minkowski spacetime or at far distance away from the
source. The first solution is identical to Eq. (16). Thus the term $m^{2}$
given by Eq. (78) is exactly same as the scalar field mass square term
$m_{\phi}^{2}$ given in Eq. (17) for our model, when $\tilde{R}=R_{0}$.
Therefore this solution suggests that there exists a massive scalar mode of
polarization of GWs in the theory apart from the massless tensor modes.

At very far distance away from the source, we can consider
$\bar{g}_{\mu\nu}=\eta_{\mu\nu}$, i.e. the Minkowski metric and the background
curvature $\tilde{R}=0$. In this case, the Ricci scalar slowly varies near
zero, i.e. $R\simeq 0+\delta R$. Hence, for the Minkowski space the Eq. (66)
can be written as

$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}.$ (80)

And to the first order of $h_{\mu\nu}$, we get

$R_{\mu\nu}=\frac{1}{2}\\!\left(\partial_{\mu}\partial_{\rho}h^{\rho}_{\nu}+\partial_{\nu}\partial_{\rho}h^{\rho}_{\mu}-\partial_{\mu}\partial_{\nu}h-\square
h_{\mu\nu}\right)\\!,$ (81) $R=\partial_{\mu}\partial_{\rho}h^{\rho\mu}-\Box
h,$ (82)

where $h=\eta^{\mu\nu}h_{\mu\nu}$. So for our model, to the first order of
perturbation, the Eq. (8) becomes

$R_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}R-\frac{(2\alpha-\pi\beta)}{\pi(\beta-1)R_{c}}\left(\partial_{\mu}\partial_{\nu}R-\eta_{\mu\nu}\,\square
R\right)=0.$ (83)

Taking the trace of this Eq. (83), we get

$(\square-m_{0}^{2})R=0,$ (84)

where

$m_{0}^{2}=\frac{\pi(\beta-1)R_{c}}{6\alpha-3\pi\beta}$

with $\alpha>0$ and $\beta>0$. This is also exactly the same mass square term
$m_{\phi}^{2}\big{|}_{R_{0}\,=\,0}$ in Minkowski space given by the Eq. (18)
for our model. Next, we introduce a variable

$\bar{h}_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}h-\frac{(2\alpha-\pi\beta)}{\pi(\beta-1)R_{c}}\,\eta_{\mu\nu}R.$
(85)

The trace of this variable is

$\bar{h}=\eta^{\mu\nu}\bar{h}_{\mu\nu}=-\,h-4\frac{(2\alpha-\pi\beta)}{\pi(\beta-1)R_{c}}R.$
(86)

Using this Eq. (86) in the variable (85) we find,

$h_{\mu\nu}=\bar{h}_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}\bar{h}-\frac{(2\alpha-\pi\beta)}{\pi(\beta-1)R_{c}}\,\eta_{\mu\nu}R.$
(87)

From Eq. (85) and Eq. (87), one can easily see that both $h_{\mu\nu}$ and
$\bar{h}_{\mu\nu}$ are interchangeable, i.e. replacing $\bar{h}_{\mu\nu}$ by
$h_{\mu\nu}$ and vice-versa in Eq. (85) gives Eq. (87). Again, under an
infinitesimal coordinate transformation, $x^{\mu}\rightarrow
x^{\mu\prime}=x^{\mu}+\varsigma^{\mu}$, we have

$h_{\mu\nu}^{\prime}=h_{\mu\nu}-\partial_{\mu}\varsigma_{\nu}-\partial_{\nu}\varsigma_{\mu}.$
(88)

The trace of this equation is

$h^{\prime}=h-2\,\partial_{\mu}\varsigma^{\mu}.$ (89)

And

$\bar{h}_{\mu\nu}^{\prime}=\bar{h}_{\mu\nu}-\partial_{\mu}\varsigma_{\nu}-\partial_{\nu}\varsigma_{\mu}+\eta_{\mu\nu}\,\partial_{\rho}\varsigma^{\rho}.$
(90)

The trace trace of this equation gives,

$\bar{h}^{\prime}=\bar{h}+2\,\partial_{\rho}\varsigma^{\rho}.$ (91)

Here we raise or lower the indices with the help of $\eta_{\mu\nu}$, i.e. with
the Minkowski metric. The Lorentz gauge condition
$\partial^{\mu}\bar{h}_{\mu\nu}^{\prime}=0$ can be obtained if
$\varsigma_{\mu}$ satisfies
$\square\varsigma_{\nu}=\partial^{\mu}\bar{h}_{\mu\nu}$. The Lorentz gauge
condition does not constrain the gauge freedom and there is always a
possibility to choose the transverse and traceless conditions, i.e.
$\partial^{\mu}\bar{h}_{\mu\nu}=0$ and
$\bar{h}=\eta^{\mu\nu}\bar{h}_{\mu\nu}=0$ Liang_2017 ; Corda_2007 ; Corda_2008
; Capozziello_2008 . By using the transverse traceless gauge condition and
substituting the Eq. (87) into Eq. (81), we get

$R_{\mu\nu}=\frac{1}{2}\\!\left[-\,\square\bar{h}_{\mu\nu}+2\,\frac{(2\alpha-\pi\beta)}{\pi(\beta-1)R_{c}}\,\partial_{\mu}\partial_{\nu}R+\frac{(2\alpha-\pi\beta)}{\pi(\beta-1)R_{c}}\,\eta_{\mu\nu}\,\square
R\right]\\!.$ (92)

Plugging Eq. (92) into Eq. (83), we obtain

$3\,\frac{2\alpha-\pi\beta}{2\pi(\beta-1)R_{c}}\,\eta_{\mu\nu}(\square-
m_{0}^{2})R-\frac{1}{2}\,\square\bar{h}_{\mu\nu}=0.$ (93)

Combining Eq.s (84) and (93), we get

$\square\bar{h}_{\mu\nu}=0,$ (94)

which is the wave equation of the massless tensor field. The solution to this
Eq. (94) is Corda_2007 ; Corda_2008

$\bar{h}_{\mu\nu}=e_{\mu\nu}\exp(iq_{\mu}x^{\mu})+\text{c.c.},$ (95)

where $\eta_{\mu\nu}q^{\mu}q^{\nu}=0$ and $q^{\mu}e_{\mu\nu}=0$. Whereas the
solution to the massive scalar field Eq. (84) is given by Corda_2007 ;
Corda_2008 ,

$R=\psi=\psi_{0}\exp(ip_{\mu}x^{\mu})+\text{c.c.},$ (96)

where $\eta_{\mu\nu}p^{\mu}p^{\nu}=-\,m^{2}_{0}$. Assuming the GWs propagation
direction along $z$, the general solution can be written as Liang_2017

$h_{\mu\nu}=\bar{h}_{\mu\nu}(t-z)+\frac{2\alpha-\pi\beta}{\pi(1-\beta)R_{c}}\,\eta_{\mu\nu}\,\psi(vt-z),$
(97)

where $\bar{h}_{\mu\nu}$ is transverse and traceless and it represents the
standard spin-2 graviton, and $\psi$ represents the scalar field which is
massive in nature and travels with a speed less than $c$. The solution of the
scalar part along $z$ axis, i.e. $\psi(vt-z)$ can be expressed as

$\psi=\psi_{0}\,e^{-i\omega t\,+\,ikz}$ (98)

and hence the mass of the field in terms of $k$ and $\omega$ is

$m_{0}=\sqrt{\omega^{2}-k^{2}}.$ (99)

###  V.2 Calculation of Exact Polarization amplitudes and Newman-Penrose
Quantities of the Model

In 1973, a powerful method was introduced in Ref. Eardley_1973  which deals
with the study of the properties of GWs in any metric theory of gravity. This
method involves analysing all the relevant components of Riemann tensor, which
results relative acceleration between two test particles. They used a null-
tetrad basis in order to calculate the Newman-Penrose quantities Newman_1962 .
In the Newman-Penrose formalism, there are ten $\Psi$’s, nine $\Phi$’s, and a
$\Lambda$, which are algebraically independent and represent the irreducible
parts of the Riemann tensor $R_{\lambda\mu\kappa\nu}$. They are known as
Newman-Penrose quantities. But in case of plane waves or nearly plane waves,
the differential and symmetry properties of $R_{\lambda\mu\kappa\nu}$ reduce
the number of independent, nonvanishing components, to six. Hence, in this
formalism, the set $\\{\Psi_{2},\Psi_{3},\Psi_{4},\Phi_{22}\\}$ is used to
describe the six independent components of GWs in the metric theory. In the
tetrad basis, the Newman-Penrose quantities of the Riemann tensor are
Eardley_1973 :

$\displaystyle\Psi_{2}$ $\displaystyle=-\,\frac{1}{6}\,R_{lklk},$ (100)
$\displaystyle\Psi_{3}$ $\displaystyle=-\,\frac{1}{2}\,R_{lkl\overline{m}},$
(101) $\displaystyle\Psi_{4}$
$\displaystyle=-\,R_{l\overline{m}l\overline{m}},$ (102)
$\displaystyle\Phi_{22}$ $\displaystyle=-\,R_{lml\overline{m}}.$ (103)

It should be noted that, $\Psi_{3}$ and $\Psi_{4}$ are complex. Therefore,
each one of them is capable of describing two independent polarizations. One
polarization mode for the real part and one for the imaginary part. Thus total
number of polarization modes is 6.

The tetrad components of Ricci tensors can be expressed as Eardley_1973 ,

$\displaystyle R_{lk}$ $\displaystyle=R_{lklk},$ (104) $\displaystyle R_{ll}$
$\displaystyle=2\,R_{lml\overline{m}},$ (105) $\displaystyle R_{lm}$
$\displaystyle=R_{lklm},$ (106) $\displaystyle R_{l\overline{m}}$
$\displaystyle=R_{lkl\overline{m}},$ (107)

and the Ricci scalar is

$R=-\,2\,R_{lk}=-\,2\,R_{lklk}.$ (108)

In normal coordinate system Eardley_1973 ,

$\displaystyle\Psi_{2}$ $\displaystyle=-\,\dfrac{1}{6}\,R_{ztzt},$
$\displaystyle\Psi_{3}$
$\displaystyle=-\,\dfrac{1}{2}\,R_{xtzt}+\dfrac{i}{2}\,R_{ytzt},$
$\displaystyle\Psi_{4}$ $\displaystyle=-\,R_{xtxt}+R_{ytyt}+2i\,R_{xtyt},$
$\displaystyle\Phi_{22}$ $\displaystyle=-\,R_{xtxt}-R_{ytyt}.$

Although the amplitudes $\\{\Psi_{2},\Psi_{3},\Psi_{4},\Phi_{22}\\}$ of a wave
depend on the observer Eardley_1973 , there are certain invariant statements
about them that hold true for all the standard observers if they hold true for
any one. These statements characterize the invariant $E(2)$ classes of waves.
For a standard observer, under the assumptions that (a) the wave travels in
the $+z$ direction, and (b) the same frequency for a monochromatic wave is
observed, the $E(2)$ classes are:

  * •

Class $II_{6}$: $\Psi_{2}\neq 0$. Standard observers measure the same non-
vanishing amplitude in the $\Psi_{2}$ mode. Presence or absence of all other
modes is observer-dependent;

  * •

Class $III_{5}$: $\Psi_{2}=0,~{}\Psi_{3}\neq 0$. Standard observers measure
the absence of $\Psi_{2}$ and the presence of $\Psi_{3}$. Presence or absence
of $\Psi_{4}$ and $\Phi_{22}$ is observer-dependent;

  * •

Class $N_{3}$: $\Psi_{2}=\Psi_{3}=0,~{}\Psi_{4}\neq 0,\Phi_{22}\neq 0$.
Presence or absence of all modes is observer-independent;

  * •

Class $N_{2}$: $\Psi_{2}=\Psi_{3}=\Phi_{22}=0;~{}\Psi_{4}\neq 0$. Observer-
independent;

  * •

Class $O_{1}$: $\Psi_{2}=\Psi_{3}=\Psi_{4}=0;~{}\Phi_{22}\neq 0$. Observer-
independent;

  * •

Class $O_{0}$: $\Psi_{2}=\Psi_{3}=\Psi_{4}=\Phi_{22}=0$. Observer-independent.
All standard observers measure no wave.

In $f(R)$ gravity, the field equation derived from the Lagrangian in metric
formalism results dynamical expressions for Ricci tensor and Ricci scalar, as
we have already seen. Using this method, the expressions for Ricci tensor and
scalar are calculated in weak field limit, i.e. far from the GWs source
considering that the GW is propagating along $z$ axis. Thus, the Ricci tensor
components corresponding to directions other than $z$ and $t$ will vanish. But
from earlier sections, we see that a massive scalar mode of polarization is
present in this model and so it is not possible to use this Newman-Penrose
formalism formalism which is developed for null waves Liang_2017 . For waves
with massive propagation mode, modified Newman-Penrose formalism has to be
applied Hyun_2019 . In modified Newman-Penrose formalism, the polarization
amplitudes as well as Newman-Penrose quantities are calculated for a massive
wave subject to proper gauge condition. Considering the monochromatic wave
solution of the form:

$\displaystyle h_{\mu\nu}=C_{\mu\nu}e^{-i\omega t+ikz}\,,$ (109)

where $\omega$ is the frequency and $k$ is the wave number. However, we have
noticed that, the modified Newman-Penrose scalars and polarization amplitudes
introduced in Hyun_2019  using Lorentz gauge condition can not distinguish
breathing mode when the scalar field becomes massless. It is due to the fact
that, they have used the transverse traceless condition to make the non-tensor
modes vanish when $\omega=k$. But $\omega=k$ can’t demand $p_{6}^{(b)}=0$
because, the breathing mode is massless in nature Gogoi_2019 . Breathing modes
satisfy the transverse condition but not the traceless condition. It implies
that a model having $\omega=k$ can have massless breathing mode of
polarization Gogoi_2019 . Keeping this fact in mind, we have modified the
polarization amplitudes and Newman-Penrose quantities. According to our
calculations, modified polarization amplitudes are expressed as:

$\displaystyle p_{1}^{(l)}\,\,$
$\displaystyle=\frac{1}{2}\\!\left(\frac{\omega^{2}-k^{2}}{\omega^{2}+k^{2}}\right)\omega^{2}(h_{tt}+h_{zz})-\frac{1}{2}\left(\omega^{2}-k^{2}\right)h_{tt}\,,$
$\displaystyle p_{2}^{(x)}\,$
$\displaystyle=\frac{1}{2}\left(\omega^{2}-k^{2}\right)h_{xz}\,,$
$\displaystyle p_{3}^{(y)}\,$
$\displaystyle=\frac{1}{2}\left(\omega^{2}-k^{2}\right)h_{yz}\,,$
$\displaystyle p_{4}^{(+)}$
$\displaystyle=\frac{1}{2}\,\omega^{2}\left(h_{xx}-h_{yy}\right)\,,$
$\displaystyle p_{5}^{(\times)}$
$\displaystyle=\frac{1}{2}\,\omega^{2}h_{xy},$ $\displaystyle p_{6}^{(b)}\,$
$\displaystyle=\frac{1}{2}\,\omega^{2}(h_{xx}+h_{yy})\,.$ (110)

Here, we have not applied the traceless condition to the breathing mode. These
are the exact polarization amplitudes of the wave. These expressions are valid
for any metric theory. The modified Newman-Penrose quantities now can be
expressed as:

$\displaystyle\Psi_{2}$
$\displaystyle=-\,\frac{1}{24}\left(\frac{\omega^{2}-k^{2}}{\omega^{2}+k^{2}}\right)\left[(3k^{2}-\omega^{2})h_{tt}+(k^{2}-3\omega^{2})h_{zz}\right]\,,$
$\displaystyle\Psi_{3}$
$\displaystyle=\frac{1}{8}\,\frac{(\omega-k)(\omega+k)^{2}}{\omega}\,(h_{xz}-ih_{yz})\,,$
$\displaystyle\Psi_{4}$
$\displaystyle=\frac{1}{8}\,(\omega+k)^{2}(h_{xx}+h_{yy})-\frac{1}{4}(\omega+k)^{2}(h_{yy}+ih_{xy})\,,$
$\displaystyle\Phi_{22}$
$\displaystyle=\frac{1}{8}\,(\omega+k)^{2}(h_{xx}+h_{yy})\,.$ (111)

These equations differ from those in Ref. Hyun_2019 .

Now, for our model the general wave solution is given by Eq. (97). Using this
Eq. (97) in the above set of Eq.s (110), we found the polarization amplitudes
for our model as

$\displaystyle p_{1}^{(l)}\,\,$
$\displaystyle=\dfrac{1}{2}\,m^{2}_{0}\,C_{1}\psi\,,$ $\displaystyle
p_{2}^{(x)}\,$ $\displaystyle=0\,,$ $\displaystyle p_{3}^{(y)}\,$
$\displaystyle=0\,,$ $\displaystyle p_{4}^{(+)}$
$\displaystyle=-\,\dfrac{1}{2}\,(\ddot{\bar{h}}_{xx}-\ddot{\bar{h}}_{yy})\,,$
$\displaystyle p_{5}^{(\times)}$
$\displaystyle=-\,\frac{1}{2}\,\ddot{\bar{h}}_{xy}\,,$ $\displaystyle
p_{6}^{(b)}\,$ $\displaystyle=\omega^{2}C_{1}\psi\,,$ (112)

where $m^{2}_{0}$ is given by Eq. (99) and
$C_{1}=\frac{2\alpha-\pi\beta}{\pi(1-\beta)R_{c}}$ (see Eq. (97)). From the
above expressions we can calculate the Newman-Penrose quantities. Note that,
above results suggest, there are $4$ non-zero polarization amplitudes in the
theory. Using Eq.s (111), we’ve calculated the Newman-Penrose quantities for
the model as

$\displaystyle\Psi_{2}$ $\displaystyle=\dfrac{1}{12}\,m^{2}_{0}\,C_{1}\psi\,,$
$\displaystyle\Psi_{3}$ $\displaystyle=0\,,$ $\displaystyle\Psi_{4}$
$\displaystyle=(\ddot{\bar{h}}_{yy}+i\ddot{\bar{h}}_{xy})\,,$
$\displaystyle\Phi_{22}$
$\displaystyle=\dfrac{1}{4}\,(\omega+k)^{2}\,C_{1}\psi\,.$ (113)

Thus the $E(2)$ classification of the model is $II_{6}$. The model exhibits
non-zero Newman-Penrose quantities for $4$ polarization modes viz., tensor
plus, tensor cross, scalar transverse massless breathing mode and scalar
longitudinal massive mode of polarization. However, the degrees of freedom
associated with the theory is $3$. This suggests that the breathing mode and
the longitudinal mode exist in a mixed state to give rise to a single
polarization mode. If $m_{0}=0$, the massive longitudinal mode will vanish,
giving $\Psi_{2}=0$. Note that $m_{0}=0$ is not a sufficient condition to
imply the absence of scalar degrees of freedom in the theory. It is because,
in $f(R)$ theory there exists massless breathing mode which is transverse but
not traceless. Absence of scalar degrees of freedom requires both $\Psi_{2}=0$
and $\Phi_{22}=0$. When both $m_{0}$ and $\psi(vt-z)$ vanish, the theory
reduces to GR giving only tensor modes of polarizations.

##  VI Detection of Polarization Modes of GWs - A Review

Experimental detection of polarization modes of GWs is very important to know
the exact nature of GWs and hence in checking the viabilities of modified
gravity theories. In this section we discuss the Pulsar Timing Arrays (PTAs)
as a tool to distinguish between different polarization modes. Moreover, we
include a discussion on the results based on our model.

PTAs play a significant role in the indirect detection of GWs. They are also
used for numerous astrophysical applications. In 1968, Counselman & Shapiro
explained that the observations of pulsars could be used to test GR
Counselman_1968 . Later in 1982, the first millisecond pulsar was discovered
Backer_1982 . Till now, a pretty good number of millisecond pulsars has been
discovered. The advantage of these pulsars over the normal pulsars is that
they are very stable. Their arrival times can be measured and predicted with a
good accuracy. This allows to use these pulsars as a probe to search for GWs.
In 2004, the Parkes Pulsar Timing Array (PPTA) project began with the Parkes
64 m telescope Manchester2013 ; Oslowski2019 . After three years, the North
American Nanohertz Observatory for Gravitational Waves (NANOGrav) in North
America was founded Brazier2019 ; Arzoumanian2016 . NANOGrav uses the Arecibo
and Green Bank telescopes to observe around 36 pulsars. In the same year, the
European Pulsar Timing Array (EPTA) project was also founded Desvignes2016 .
Using the Sardinian, Effelsberg, Nancay, Westerbork and Jodrell Bank
telescopes, EPTA observes around 42 pulsars. Later, by combining these three
projects, the International Pulsar Timing Array (IPTA) was formed Verbiest2016
; Perera2019 . In this study, we have used some selected data from PPTA
Manchester2013  and IPTA Verbiest2016 ; Perera2019 .

Figure 16: Variation of correlation functions $(\zeta(\theta))$ with respect
to $\theta$. Plot on the left shows the correlation function for tensor modes
and massless breathing mode of polarization of GWs. Plot on the right shows
the correlation function for the longitudinal mode of polarization for
different values of the mass of the mode as predicted by the model (1) along
with correlation functions for some selected pulsars obtained from PPTA
Manchester2013  and IPTA Verbiest2016 ; Perera2019  data.

From the Ref. Lee_2008 , we can have the correlation functions for different
polarization modes. The calculation of correlation functions for the tensor
and breathing modes are model independent, but in the case of massive
longitudinal mode, the correlation function is model dependent as it depends
on the mass of the scaler graviton. The correlation function for tensor modes
is Lee_2008

$C^{+,\times}(\theta)=\xi^{GR}(\theta)\int_{0}^{\infty}\dfrac{|h_{c}^{+,\times}|^{2}}{24\pi^{2}f^{3}}\,df,$
(114)

where

$\xi^{GR}(\theta)=\dfrac{3\,(1-\cos\theta)}{4}\log\left(\dfrac{1-\cos\theta}{2}\right)+\dfrac{1}{2}-\dfrac{1-\cos\theta}{8}+\dfrac{\delta(\theta)}{2},$

and $\theta$ is the angular separation between two pulsars. For the scalar
modes it is Lee_2008

$C^{b}(\theta)=\xi^{b}(\theta)\int_{0}^{\infty}\dfrac{|h_{c}^{b}|^{2}}{12\pi^{2}f^{3}}\,df,$
(115)

where

$\xi^{b}(\theta)=\dfrac{1}{8}\big{[}\cos\theta+3+4\,\delta(\theta)\big{]}.$

The normalized correlation function in general is given by:

$\zeta(\theta)=\dfrac{C(\theta)}{C(0)}.$

These are the correlation functions for tensor modes and massless breathing
mode of GW polarization. But in our model, there exists a massive longitudinal
mode. Thus to see the effect of massive mode, we follow the Refs. Jenet_2005 ;
Lee_2013 ; Lee_2010  in which the timing residual induced by GWs is expressed
as

$R=-\,\frac{1}{\cal S}\,A^{ij}H_{ij},$ (116)

here ${\cal S}=2\left(1+({c}/{{\omega_{g}}}){{\mathbf{k}_{g}}}\cdot\hat{\rm\bf
n}\right)$ gives the dispersion relation of the GWs, the terms $\omega_{g}$
and $k_{g}$ connect the mass of the longitudinal mode of polarization $m_{g}$
by the relation $m_{g}^{2}=\omega_{g}^{2}-k_{g}^{2}$, $A^{ij}\equiv\hat{\rm\bf
n}^{i}\hat{\rm\bf n}^{j}$ and
$H_{ij}=\int_{0}^{\tau}h_{ij}(\tau,0)-h_{ij}(\tau-|{\rm\bf D}|/c,{\rm\bf
D})\,d\tau$. $\rm\bf D$ is the displacement vector from the observer to the
pulsar and $\hat{\rm\bf n}^{i}$ and $\hat{\rm\bf n}^{j}$ are two unit vectors
pointing to two pulsars. The correlation coefficient $C$ between two different
pulsars is given by

$C_{1,2}(\theta)=\langle R_{1}R_{2}\rangle={A}_{1}{A}_{2}\langle{\cal
S}_{1}{\cal S}_{2}{H}_{1}{H}_{2}\rangle,$

where the sub-scripts are indices for the pulsars. With these assumptions and
following Ref. Lee_2008 , the correlation functions are calculated numerically
for different values of $m_{g}$ (see Fig. 16). In terms of our model,

$m_{g}^{2}\equiv
m^{2}_{\phi}=\left[\frac{R_{c}e^{R/R_{c}}\left(\pi\left(R^{4}+R_{c}^{4}\right)^{2}-8\alpha
R^{5}R_{c}^{3}\right)-\pi\beta(R+R_{c})\left(R^{4}+R_{c}^{4}\right)^{2}}{3\pi\beta\left(R^{4}+R_{c}^{4}\right)^{2}-6\alpha
R_{c}^{4}e^{R/R_{c}}\left(R_{c}^{4}-3R^{4}\right)}\right]_{R\,=\,R_{0}}\\!\\!\\!\\!\\!.$

All these numerically calculated correlation functions are plotted with
respect to $\theta$ in Fig. 16. In this figure we have used three values of
$m_{g}$ for the following sets of parameters:

(i) $m_{g}=0$ for $R=0~{}m^{-2},R_{c}=0~{}m^{-2},~{}\alpha,\beta\in[0,1]$,

(ii) $m_{g}=10^{-24}$ eV for
$R=0~{}m^{-2},R_{c}=10^{-50}~{}m^{-2},\alpha=0.05,\beta=0.031831$ and

(iii) $m_{g}=10^{-25}$ eV for $R=4.44\times 10^{-52}~{}m^{-2},R_{c}=5\times
10^{-51}~{}m^{-2},\alpha=0.0314565,\beta=0.17268558.$

To check the experimental viability of our model, we have calculated the
correlation functions for GWs using some selected data from IPTA and PPTA data
set Manchester2013 ; Verbiest2016 ; Perera2019 ; Hobbs2006 ; Edwards2006  as
mentioned above, which are also plotted in Fig. 16. Although we do not have a
clear conclusion, which requires more observation period as well as data, we
can still see that these experimental data could not directly rule out the
existence of extra polarization modes of GWs. However, to distinguish between
polarization modes, we need to wait for more PTA data with GW events.

##  VII Summary and Conclusion

In this work, we have introduced a new $f(R)$ gravity toy model and studied
the polarization modes of GWs in it. The study shows that, in metric
formalism, there exists $3$ polarization modes of GWs viz., tensor plus mode,
tensor cross mode and scalar mode. The scalar mode is a mixed state of
massless breathing mode and massive longitudinal mode. The tensor modes of
polarization are transverse, traceless and massless in nature. The scalar
breathing mode is transverse but exists with non vanishing trace and massless
in nature. On the other hand, the scalar longitudinal mode is massive in
nature and hence propagates with speed less than that of tensor modes. When
the scalar field becomes massless, the longitudinal mode vanishes and only the
massless scalar breathing mode exists in the scalar degrees of freedom. As an
experimental correspondence of our model prediction on the modes of GWs, we
have compared the correlation function of massive longitudinal mode as
predicted by the model with that of the some selected PTAs data of PPTA and
IPTA. The result is found to be quite encouraging.

We have also shown that a wisely selected set of the parameters easily allows
the model to pass the solar system tests, which is a very important
requirement for the viability of a model. Further, the model has been
constrained using combined CMB, BAO and $\sigma_{8}-\Omega_{m}$ relationship
from the PSZ catalog and Abell 1689 galaxy cluster data. The model can be
easily constrained with the help of the parameters $R_{c}$, $\alpha$ and
$\beta$. Further, we have also constrained the model using the GW event
GW170817. It is seen that, the model can withstand the constraints put by
GW170817 and hence can be included as a post GW170817 viable model.

In the present work, we have studied only few properties of the model.
However, for a proper understanding of the model characteristics, a detailed
study is required. As such, in future the model can be checked also for
various stabilities and constraints as well as in different cosmological and
astrophysical contexts, which will give us more information about the
viability of the model. Moreover, it is to be noted that in the Palatini
formalism the polarization modes of GWs in $f(R)$ gravity is model
independent. So, this formalism can not distinguish our model from other
$f(R)$ gravity models in this context. Nevertheless, there may be some
cosmological variations and stellar structure differences of the model in
Palatini formalism, which might be useful to study the generation of GWs in
such situations.

## Acknowledgement

A part of this work was done during a visit of authors to IUCAA, Pune. Authors
are grateful to IUCAA for the hospitality during their stay.

## References

  * (1) A. G. Riess et al., Astron. J. 116, 1009 (1998).
  * (2) S. Perlmutter et al., Astrophys. J. 517, 565 (1999).
  * (3) N. A. Bahcall, J. P. Ostriker, S. Perlmutter and P. J. Steinhardt, Science 284, 1481 (1999) [arXiv:astro-ph/9906463].
  * (4) R. P. Kirshner, PNAS 96, 4224 (1999).
  * (5) C. Corda, Int. J. Mod. Phys. A 23 (10), 1521 (2008).
  * (6) B. P. Abbott et al., Phys. Rev. Lett. 116, 061102 (2016).
  * (7) A. A. Starobinsky, JETP Lett. 86, 157 (2007) [arXiv:0706.2041].
  * (8) W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007) [arXiv:0705.1158].
  * (9) S. Tsujikawa, Phys. Rev. D 77, 023507 (2008) [arXiv:0709.1391].
  * (10) P. Zhang, Phys. Rev. D 73, 123504 (2006) [arXiv:astro-ph/0511218].
  * (11) G. Cognola et al., Phys. Rev. D 77, 046009 (2008) [arXiv:0712.4017].
  * (12) E.V. Linder, Phys. Rev. D 80, 123528 (2009) [arXiv:0905.2962].
  * (13) C. Corda, JCAP 04, 009 (2007).
  * (14) D. Liang, Y. Gong, S. Hou and Y. Liu, Phys. Rev. D 95, 104034 (2017).
  * (15) D. J. Gogoi and U. D. Goswami, arXiv:1901.11277 (2019).
  * (16) Y.-C. Chen, C.-Q. Geng, C.-C. Lee, and H. Yu, Eur. Phys. J. C 79, 93 (2019).
  * (17) L. Amendola, S. Tsujikawa, Dark Energy: Theory and Observations (Cambridge University Press, Cambridge, 2015).
  * (18) L. Amendola, R. Gannouji, D. Polarski, and S. Tsujikawa, Phys. Rev. D 75, 083504 (2007) [arXiv:gr-qc/0612180].
  * (19) K. Bamba, C.Q. Geng and C.C. Lee, JCAP 1011, 001 (2010) [arXiv:1007.0482].
  * (20) H. Motohashi, A. A. Starobinsky, and J. Yokoyama, JCAP 06, 006 (2011) [arXiv:1101.0744].
  * (21) S. Capozziello, A. Stabile and A. Troisi, Class. Quantum Grav. 24, 2153 (2007) [arXiv:gr-qc/0703067].
  * (22) P. Teyssandier and P. Tourrenc, J. Math. Phys. (N.Y.) 24, 2793 (1983).
  * (23) U. D. Goswami and K. Deka, Int. J. Mod. Phys. D 22, 1350083 (2013) [arXiv:1303.5868].
  * (24) S. Chakraborty, S. Pal, and A. Saa, Phys. Rev. D 99, 024020 (2019) [arXiv:1812.01694].
  * (25) J.-Q. Guo, Int. J. Mod. Phys. D 23, 1450036 (2014) [arXiv:1306.1853].
  * (26) L. Boubekeur et al., Phys. Rev. D 90, 103512 (2014).
  * (27) S. Desai, Phys. Lett. B 778, 325 (2018) [arXiv:1708.06502].
  * (28) S. Gupta and S. Desai, Class. Quantum Grav. 36, 105001 (2019) [arXiv:1811.09378].
  * (29) S. Jana and S. Mohanty, Phys. Rev. D 99, 044056 (2019) [arXiv:1807.04060].
  * (30) T. Chiba, T. L. Smith, and A. L. Erickcek, Phys. Rev. D 75, 124014 (2007).
  * (31) M. Cataneo et al., Phys. Rev. D 92, 044009 (2015) [arXiv:1412.0133].
  * (32) R. C. Nunes, S. Pan, E. N. Saridakis, and E. M. C. Abreu, J. Cosmol. Astropart. Phys. 01, 005 (2017).
  * (33) P. A. R. Ade et al. (Planck Collaboration), A$\&$A 571, A20 (2014) [arXiv:1303.5080].
  * (34) C. Corda, J. Cosmol. Astropart. Phys. 04, 009 (2007).
  * (35) C. Corda, Int. J. Mod. Phys. A23, 1521 (2008).
  * (36) S. Capozziello, C. Corda, and M. F. De Laurentis, Phys. Lett. B 669, 255 (2008).
  * (37) D. M. Eardley, D.L. Lee, A.P. Lightman, Phys. Rev. D 8, 3308 (1973).
  * (38) E. Newman, R. Penrose, J. Math. Phys. 3, 566 (1962).
  * (39) Y.-H. Hyun, Y. Kim, and S. Lee, Phys. Rev. D 99, 124002 (2019).
  * (40) C. C. Counselman III and Shapiro II, Science 162, 352 (1968).
  * (41) D. C. Backer, S. R. Kulkarni and C. Heiles et al., Nature 300, 615 (1982).
  * (42) R. N. Manchester et al., Publ. Astron. Soc. Aust. 30, e017 (2013) [arXiv:1210.6130].
  * (43) S. Osłowski et al., Mon. Not. R. Astron. Soc. 488, 868 (2019) [arXiv:1906.09793].
  * (44) A. Brazier et al., arXiv:1908.05356 (2019).
  * (45) Z. Arzoumanian et al., ApJ 821, 13 (2016) [arXiv:1508.03024].
  * (46) G. Desvignes et al., Mon. Not. R. Astron. Soc. 458, 3341 (2016) [arXiv:1602.08511].
  * (47) J. P. W. Verbiest et al., Mon. Not. R. Astron. Soc. 458, 1267 (2016) [arXiv:1602.03640].
  * (48) B. B. P. Perera et al., Mon. Not. R. Astron. Soc. 490, 4666 (2019) [arXiv:1909.04534].
  * (49) K. J. Lee, F. A. Jenet, and R. H. Price, ApJ 685, 1304 (2008).
  * (50) F. A. Jenet, G. B. Hobbs, K. J. Lee, and R. N. Manchester, ApJ 625, L123 (2005).
  * (51) K. J. Lee, Class. Quantum Grav. 30, 224016 (2013) [arXiv:1404.2090].
  * (52) K. Lee, F. A. Jenet, R. H. Price, N. Wex, and M. Kramer, ApJ 722, 1589 (2010) [arXiv:1008.2561].
  * (53) G. B. Hobbs, R. T. Edwards, and R. N. Manchester, Mon. Not. R. Astron. Soc. 369, 655 (2006) [arXiv:astro-ph/0603381].
  * (54) R. T. Edwards, G. B. Hobbs, and R. N. Manchester, Mon. Not. R. Astron. Soc. 372, 1549 (2006) [arXiv:astro-ph/0607664].