---
abstract: 'We’ve investigated the different polarization modes of Gravitational Waves (GWs) in $f(R)$ gravity power law model in de Sitter space. It is seen that the massive scalar field polarization mode exists in this model. The mass of the scalar field depends highly on the background curvature and the power term $n$. However, we found that the model doesnot exhibit a massive scalar mode for $n=2$ and instead it shows a breathing mode in addition to the tensor plus and cross modes. Thus mass of the scalar field is found to vary with $n$ within the range $1 < n < 2$.'
author:
- Dhruba Jyoti Gogoi$^1$ and Umananda Dev Goswami$^1$
title: 'Gravitational Waves in $\mathbf{f(R)}$ Gravity Power Law Model'
---

Introduction
============

Inspired by scientific curiosity, modification of Einstein’s gravity was at first attempted by Weyl in 1919 [@weyl_neue_1919] and thereafter by Eddington in 1923 [@eddington_mathematical_1924] by including higher order curvature invariants in gravitational action. It is seen that these modifications had no experimental motivations and they were almost ignored up to 1962, when it was for the first time realized that modification of gravitational action had some merits. Gravity from General Relativity (GR) is not renormalizable and thus it is not possible to quantize it according to conventional methods. In 1962, Utiyama and DeWitt showed that renormalization at one loop demands that the Einstein-Hilbert action be modified by higher order curvature terms [@utiyama_renormalization_1962]. Later, a number of drawbacks of GR have been observed including the inability to explain cosmic acceleration. To overcome these drawbacks of GR a number of new theories have been proposed suggesting the further modification of GR. Among them, $f(R)$ gravity is a prominent one keeping itself standing against the challenges and problems till now [@sotiriou_fr_2010]. Another important reason to choose $f(R)$ gravity is its less complexity compared to the other modified theories of gravity.

After around 100 years of Einstein’s prediction of Gravitational Waves (GWs), the Laser Interferometer Gravitational Wave Observatory (LIGO) Scientific Collaboration and Virgo Collaboration announced the detection of GWs for the first time in September 14, 2015, which was named as GW150914 [@abbott_observation_2016]. This observation was later followed by four other events, viz. GW151226 [@abbott_gw151226_2016], GW170104 [@abbott_gw170104_2017], GW170814 [@abbott_gw170814:_2017] and GW170817 [@abbott_gw170817:_2017]. These events opened new directions in testing GR and modified theories of gravity. Therefore the study of GWs in modified gravity will play a very important role in modifying GR and on the predictions of GR. The experimental results obtained from the sector of GWs can effectively constrain the modified gravity models to a fair extent and can even contribute to the bottom-up approach in forming or proposing modified gravity theories and models in an efficient way.

In the metric formalism of $f(R)$ gravity, the GWs have other modes of polarizations besides the tensor modes of polarization found in GR. The propagating degrees of freedom of GWs in $f(R)$ gravity is found to be 3 [@myung_propagating_2016] by Myung and later it was confirmed in D. Liang et. al. [@liang_polarizations_2017]. D. Liang et. al. in their extensive study on the Starobinsky model showed that there exists a mixed state of massless breathing mode and massive scalar mode of polarization besides the tensor modes. Thus the Newmann Penrose (NP) formalism cannot be applied in that study due to the existence of massive scalar modes of GW polarization. In our study, we’ll restrict ourselves to the power law model of $f(R)$ gravity given by $f(R) = \alpha R^n$, where $\alpha$ is an arbitrary constant and $n$ is a real number. This model has been pursued in [@jaime_about_2013] and where it was concluded that the model can produce a late time acceleration in the Universe or a transient matter domination era, but not both. Same results were also obtained in the Refs. [@amendola_power-laws_2007; @amendola_are_2007; @amendola_conditions_2007] previously. In the Refs. [@capozziello_quintessence_2003; @carloni_cosmological_2005; @goswami_cosmological_2013] this model was discussed for the different special cases along with $n=-1$ and $n=3/2$, in which $n=3/2$ case is conformally equivalent to Liouville field theory. In another paper [@faraoni_r_2011], the model was used for invoking the chameleon mechanism to work within the solar system and in which it was shown that for this to happen the value of $n$ is required to be very close to $1$. The result rendered from the work implied that the model is a poor candidate for a realistic alternative to dark energy. The cosmological dynamics of the model also has been studied in Ref. [@carloni_cosmological_2005].

In this work, we’d like to test the model in the GW regime to find out the possibility of its viability and the drawbacks. Another important point of this model is that it is the only model in $f(R)$ theories in the metric formalism which can give pure massless breathing mode as a special case. Here we’ll study the model with the special case of $n=2$ along with the stability and polarization modes.

The paper is organized as follows. In the next section the variation of mass of the scalar field in de Sitter space with the power term $n$ is studied. In the third section stability of the model is studied in de Sitter space. In the fourth section, the contents of the theory and the model are studied with the help of operator formalism. In there, we tried to show the contents of the model for a specific case of $n=2$ and we show that in $f(R)$ theory of modified gravity a pure massless scalar field is also possible as in case of Scalar Tensor theory. The fifth section contains a study towards the equivalence of the theory with the Scalar-Tensor theory and the required conditions for such equivalence. The tensor and the scalar polarization modes are studied with the help of geodesic deviation in the sixth section and with the help of Newman-Penrose formalism in the seventh section. In section eight, we have discussed a possible way to check the validity of the model experimentally. 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.

Scalar and tensor fields from the model {#sec2}
=======================================

In this study of GWs we use vacuum space. The action of $f(R)$ gravity in vacuum is given by $$S =\dfrac{1}{2\kappa} \int d^4x \sqrt{-g}\, f(R),
\label{eq1}$$ where $f(R)$ is a function of Ricci curvature $R$. The vacuum field equations obtained from this action are given by $$\label{eq2}
 f'(R)R_{\mu\nu} -\dfrac{1}{2} f(R)g_{\mu\nu} - \nabla_\mu \nabla_\nu f'(R) 
 + g_{\mu\nu}\, \square f'(R) = 0,$$ where $\square = g^{\mu\nu}\,\nabla_\mu \nabla_\nu$ and the prime over $f(R)$ denotes the derivative with respect to $R$. Taking trace of vacuum field equation (\[eq2\]), we obtain $$\label{eq3}
f'(R) R +3\, \square f'(R) - 2 f(R) = 0.$$

Now, we assume that there is a propagating GWs in the spacetime, which will perturb the metric around its background value. Since the perturbation is usually very small and hence if we consider the background metric as $\bar{g}_{\mu\nu}$ then to the first order of perturbation value $h_{\mu\nu}$, we may express the spacetime metric as $$g_{\mu\nu}=\bar{g}_{\mu\nu} + h_{\mu\nu},\;\; \mbox{where}\; \vert h_{\mu\nu} \vert << \vert \bar{g}_{\mu\nu} \vert.
\label{eq4}$$ In view of this perturbation, expanding the Ricci tensor and the Ricci scalar upto the first order of $h_{\mu\nu}$, we may write $$\begin{aligned}
\notag
 R_{\mu\nu}& \simeq \bar{R}_{\mu\nu} + \delta R_{\mu\nu} + O(h^2)\\ 
&= \bar{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}) + O(h^2)
\label{eq5}\end{aligned}$$ and $$\begin{aligned}
\notag
 R &\simeq \bar{R} + \delta R + O(h^2)\\
 & =\bar{R} - \square h + \nabla^\mu \nabla^\nu h_{\mu\nu} -\bar{R}_{\mu\nu} h^{\mu\nu} + O(h^2),
\label{eq6}\end{aligned}$$ where $\bar{R}$ is the background curvature. Thus, due to this perturbation the trace equation (\[eq3\]) can be rewritten as $$3 f''(\bar{R}) \square \delta R + \left[f''(\bar{R})\bar{R} +3\square f''(\bar{R})- f'(\bar{R})\right] \delta R = 0.
\label{eq7}$$ Using our model $$f(R)=\alpha R^n %, $f'(R)=\alpha n R^{n-1} $, and $f''(R)=\alpha n(n-1) 
%R^{n-2}$ 
\label{eq8}$$ the equation (\[eq7\]) can be expressed in the following compact form: $$\label{eq9}
 (\square - m^2) \delta R = 0,$$ where $$\label{eq10} 
m^2 = \dfrac{(2-n)}{3(n-1)} \bar{R}.$$ The equation (\[eq9\]) is the conventional form of the scalar field equation (Klein-Gordon) with $\delta R$ as the scalar field and hence here, $m$ can be identified as the mass of this scalar field. Moreover, since this scalar field corresponding to $\delta R$ is due to ripple in spacetime, so the field can be considered as the quantized massive gravitational field.

From the expression (\[eq10\]) we see that the mass term of the scalar field depends highly on the exponent term $n$ and the background curvature $\bar{R}$. This indicates that near the massive objects like neutron stars etc., where the background curvatures are very large, the mass of the gravitational field, i.e. graviton will increase and consequently this will slow down the propagation speed of the massive mode of GWs in such regions according to this particular model. Further, for a massive graviton or massive mode of GWs $n$ should lie in between $1$ and $2$ beyond which it will correspond to a tachyonic graviton. When $n=1$ the model will give the Einstein’s case. When $n$ starts to increase from $1$ to $2$ mass of the field decreases monotonically and finally, for $n=2$ the massive scalar mode of polarization will vanish. i.e. for the value $n=2$ model will give a massless graviton. The variation of the scalar field mass with respect to $n$ for the unit value of the background curvature is shown in the Fig.\[fig1\]. One interesting observation from the mass term is that when $n$ approaches to $1$ the mass of the scalar field increases very rapidly and for $n=1$ the mass becomes infinite and hence debarred from propagating in the spacetime.

![Variation of mass $m$ of the scalar field with respect to $n$ for the unit value of the background curvature $\bar{R}$.[]{data-label="fig1"}](mass_g.eps)

The tensor field equation (\[eq2\]) in vacuum for the power law model (\[eq8\]) takes the from: $$n R^{n-1} R_{\mu\nu} -\dfrac{1}{2} g_{\mu\nu}  R^n - n \nabla_\mu \nabla_\nu R^{n-1} + n\,g_{\mu\nu} \square R^{n-1}  = 0.
\label{eq11}$$ Perturbing this field equation around the de Sitter curvature $\bar{R}$ (see the equation (\[eq5\])) we may write, $$\begin{aligned}
\nonumber
n\bar{R}\bar{R}_{\mu\nu} +n(n-1)\bar{R}_{\mu\nu} \delta R -\dfrac{1}{2} \bar{g}
  _{\mu\nu} \bar{R}^2 +n\bar{R}\,\delta R_{\mu\nu}-\dfrac{1}{2} n\, \bar{g}_{\mu\nu}\bar{R}\,\delta R\\
-\,n(n-1)\nabla_\mu \nabla_\nu(\delta R)+n(n-1)\bar{g}_{\mu\nu}\square 
  \delta R = 0.
\label{eq12}\end{aligned}$$ In the ideal case, i.e. when $R=\bar{R}$, the de Sitter curvature should satisfy the above field equation. Thus for the ideal case the above equation gives, $$\label{eq13}
  n\bar{R}_{\mu\nu}=\dfrac{1}{2}\bar{g}_{\mu\nu}\bar{R}.$$ Using this relation (\[eq13\]), we may rewrite the field equation (\[eq12\]) as $$\bar{R}\delta R_{\mu\nu} - \bar{R}_{\mu\nu} \delta R - (n-1)\nabla_\mu \nabla_\nu\delta R + (n-1)\,\bar{g}_{\mu\nu}\,\square\delta R = 0.
 \label{eq14}$$ Trace of this equation is $$\label{eq15}
  \square\delta R = 0.$$ This is the equation of massless scalar field for the case of $n=2$ (refer to equation (\[eq9\])). This is an interesting result, which shows that if the de Sitter curvature has to satisfy the field equation (\[eq12\]) then the model (\[eq8\]) has to take the value $n=2$, consequently, the massive scalar mode of polarization vanishes and a pure massless scalar mode of polarization (also known as breathing mode of polarization) is obtained.

Thus the equation (\[eq15\]) can be treated as a stability point for the model in de Sitter space which identically satisfies the stability condition discussed in the following section without modifying the Hubble constant in de Sitter space.

Again, since $n=2$ in the model (\[eq8\]) (we’ll refer this model with $n=2$ as pure $R^2$ model) makes the scalar field massless and independent of the the background curvature $\bar{R}$, so for simplicity we take this advantage to choose $\bar{R} = 0$ and perturb the tensor field equation for pure $R^2$ model (see the equation (\[eq2\])) to give $$\label{eq16}
R\delta R_{\mu \nu} - \dfrac{1}{4} R^2 \bar{g}_{\mu\nu} - \partial_\mu \partial_\nu R+\bar{g}_{\mu \nu}\square R = 0.$$ Now we define a parameter, $$\label{eq17}
h_{\mu \nu} = \bar{h}_{\mu\nu}-\dfrac{1}{2}\bar{g}_{\mu\nu}\bar{h}-\bar{g}_{\mu\nu}\delta R,$$ where taking the trace we find $\bar{h}$ as $$\bar{h} = -\, h - 4\delta R,$$ and substituting $\bar{h}$ in the definition (\[eq17\]), $\bar{h}_{\mu \nu}$ can be expressed as $$\bar{h}_{\mu \nu} = h_{\mu\nu}-\dfrac{1}{2}\bar{g}_{\mu\nu}h-\bar{g}_{\mu\nu}\delta R.$$ Choosing transverse traceless gauge condition: $ \partial^\mu \bar{h}_{\mu\nu}=0$, $\bar{h}=0$ and using the definition (\[eq17\]) in equation (\[eq5\]), we get $$\label{eq20}
\delta R_{\mu\nu} = \dfrac{1}{2} \left[-\,\square\bar{h}_{\mu\nu} + 2 \partial_\mu \partial_\nu (\delta R)
+ \bar{g}_{\mu\nu} \square \delta R  \right].$$ Finally, using this equation (\[eq20\]) and equation (\[eq15\]) in the equation (\[eq16\]), we may write $$\label{eq21}
\square \bar{h}_{\mu\nu} = 0,$$ which is the usual equation that gives the tensor modes, i.e. plus and cross modes of polarization of GWs. The solution of this equation is $$\label{eq22}
\bar{h}_{\mu\nu} = e_{\mu\nu} \exp(iq_\mu x^\mu) + c.c.,$$ where $\bar{g}_{\mu\nu} q^\mu q^\mu = 0$ and $q^\mu e_{\mu\nu}=0.$ Thus, this clearly shows that the pure $R^2$ gravity model gives three polarization modes of GWs, all massless, two of them are tensor modes while the other is the scalar or the pure breathing mode as disccussed above.

Stability of the Model
======================

According to John D. Barrow and A. C. Ottewill [@barrow_stability_1983] the condition for the existence of de Sitter solution for a constant curvature $\bar{R}$ is $$\label{eq23}
\bar{R} f'(\bar{R}) = 2 f(\bar{R}).$$ In Einstein’s case $\bar{R}=4\Lambda$ [@barrow_stability_1983], where $\Lambda$ is the cosmological constant. For our model (\[eq8\]), we get from the above condition that $$\label{eq24}
n \bar{R}^n = 2 \bar{R}^n.$$ This condition is similar to the condition obtained in equation (\[eq13\]), giving $n=2$ and which transforms the equation (\[eq9\]) into the equation (\[eq15\]) with mass for the scalar field $m=0$ as mentioned in the previous section. However, there is another possible solution for this equation (\[eq17\]) is $\bar{R}= 0$, which will lead to Minkowski spacetime having zero background curvature.

In [@faraoni_stability_2005] it was shown that the model (\[eq8\]) produce power law inflation $ a\propto t^{\bar{\alpha}} $, with $$\bar{\alpha} = \dfrac{(n-1)(1-2n)}{(n-2)}$$ for generic values of $n\neq 0,1/2,1$ and as we have already seen that this model does not allow or admit de-Sitter solution for $n\neq2$. This will be a case only when a cosmological constant is added to the action corresponding to a term with $n = 0$ [@barrow_stability_1983; @faraoni_stability_2005]. But for any positive value of $n$ Minkowski space is a solution without any cosmological constant [@faraoni_stability_2005]. The stability condition for this model from [@faraoni_stability_2005] reads $$\label{25}
\bar{R}\dfrac{2-n}{n(n-1)}\geq 0.$$ This condition gives, $1<n\leq 2$ and $n< 0$. For $n=2$ it is identically satisfied without imposing any constraints on the Hubble parameter in de Sitter space and also shows that for any values of $n$ the Minkowski space is stable [@faraoni_stability_2005]. Thus our previous results are supported by these inferences.

Contents of the Model
=====================

One of the very usefull ways to see the contents of the model is by using operator formalism. In this formalism the most general Lorentz and parity invariant Lagrangian for a tensor field $h_{\mu\nu}$ can be expressed in terms of space operators. In the de Sitter spacetime, by inverting the space operators acting on the term $h_{\mu\nu}$ in the linearized equation of motion, one can find out the propagators [@nunez_content_2004]. Thus in particular, using spin projection operators, it is possible to write the propagators in the following form:

$$\begin{gathered}
\label{eq26}
-\dfrac{P^2}{(\square + R/2)f'(R)} - \dfrac{P^1_m}{(R/2)f'(R)} - \dfrac{f'(R)/2 + R/
4 f''(R)}{(\square +R/2)f'(R)[f''(R)(3\square + R)-f'(R)]} P^0_s \\
- \dfrac{\sqrt{3}\big[ f'(R)/2 - (\square + R/4)f''(R) \big]}{(\square +R/
2)f'(R)\big[ f''(R)(3\square + R) - f'(R) \big]}(P^{0}_{sw}+P^{0}_{ws})\\
+\dfrac{\big[(4\square+R/2)f'(R)-3(2\square+R/2)^2f''(R)\big]}{(\square+R/
2)Rf'(R)\lbrace f''(R) (3\square+R)-f'(R) \rbrace} P^{0}_{w}.\end{gathered}$$

Here $P^2, P^1_m, P^0_s, P^0_{sw}, P^0_{ws}, P^0_w $ are the spin projection operators (See the Appendix and Refs. [@van_nieuwenhuizen_ghost-free_1973; @van_nieuwenhuizen_supergravity_1981] for details). As a special case, for the Minkowski spacetime ($R=0$) using $f(R)=R$ in this equation (\[eq26\]), we can have the Einstein’s case which is a standard result for the graviton propagator, given by $$\label{eq27}
-\dfrac{P^2}{\square}+ \dfrac{1}{2}\dfrac{P^0_s}{\square},$$ where $ P^0_s $ is a ghost term which cancels the mass degrees of freedom of the massive spin 2 projector $ P^2$.

To study the particle contents of the model, the terms containing $ P^0_s $ and $ P^2 $ are required to observe. From (\[eq26\]) these two terms are $$\label{eq28}
  -\dfrac{P^2 - \dfrac{1}{2}P^0_s}{(\square + R/2)f'(R)}  - \dfrac{P^0_s}
  {(\square +R/3 -f'(R)/3f''(R))f'(R)}.$$ The first term of this propagator (\[eq28\]) gives the Einstein’s case just discussed above. The second term of it for the pure $R^2$ model takes the following form: $$\label{eq29}
- \dfrac{P^0_s}{2\, \alpha\, \square R },$$ which shows the existence of massless scalar field in the model. Similarly, in general for the model (\[eq8\]) this term gives $$- \dfrac{P^0_s}{\,n\,\alpha \,(\square - m^2)\, R^{\,n-1}},$$ where $m$ is the mass term same as equation (\[eq10\]). These results justify and validate our earlier results found in the section \[sec2\]. Moreover, to ensure whether the case $n=2$ results a pure massless breathing mode of polarization, we’d study the polarization modes of the model in details in the following sections.

Equivalence with the Scalar Tensor Theory
=========================================

A straight forward way to study $f(R)$ gravity is to see its equivalence with the Scalar Tensor Theory (STT) [@chiba_1/r_2003]. This is reasonable because as like $f(R)$ theories, STTs also have two types of polarization modes, tensor modes and scalar mode. By introducing a scalar field $\phi \equiv R$, the general form of the action (\[eq1\]) can be given as $$\label{eq30}
S=\dfrac{1}{2\kappa}\int d^4 x \sqrt{-g}\,\big[\psi(\phi)R-V(\phi)\big],$$ where $\psi (\phi) = f'(\phi)$ and $v(\phi) = \phi f'(\phi)- f(\phi).$ Now, differentiating equation (\[eq30\]) with respect to $\phi$, we find $$f''(R)(R-\phi)=0.$$ This shows that when $f''(R)\neq 0$, we must have $R=\phi$. Thus $f''(R)\neq 0$ is a very important condition required to be satisfied in order to compare the $f(R)$ theory with the STT. The pure $R^2$ model satisfies this condition. In metric formalism, the STT with Brans-Dicke parameter $\omega =0$ is equivalent to the $f(R)$ gravity.

Again the action (\[eq30\]) can be written as $$\label{31}
 S=\dfrac{1}{2\kappa}\int d^4 x \sqrt{-g}\big[f(\phi)+f'(\phi)(R-\phi)\big].$$ The field equation obtained from the above action is given as $$\label{eq32}
 G_{\mu\nu}=\dfrac{1}{f'(\phi)}\big[\nabla_\mu\nabla_\nu f'(\phi)-\bar{g}_{\mu\nu}\,\square 
 f'(\phi) + \dfrac{1}{2}\,\bar{g}_{\mu\nu}\lbrace f(\phi)-\phi f'(\phi)\rbrace\big].$$ Trace of the above equation is $$\square f' = \dfrac{2}{3} f(\phi) - \dfrac{1}{3} \phi f'(\phi).$$ For the pure $R^2$ model the equation (\[eq32\]) gives $$\label{eq34}
 G_{\mu\nu}=\phi^{\,-1}\big(\partial_\mu\partial_\nu \phi-\bar{g}_{\mu\nu}\,\square \,\phi -\dfrac{1}{4} \bar{g}_{\mu\nu}\,\phi^2\big).$$ And from the trace of this equation we get, $$\label{eq35}
 \square\,\phi=0.$$ This equation is similar to the equation (\[eq15\]) which gives a massless breathing mode of polarization. In Minkowski space $\square\,\phi \equiv \square\,\delta\phi$, since $\square\,\bar{\phi} = 0$ in this space, so the solution of this equation can be written as, $$\label{eq36}
 \delta\phi=\phi_0 \exp(ip_\mu x^\mu) + c.c.,$$ where $\bar{g}_{\mu\nu}p^\mu p^\nu=0.$ Now, combining the solutions for the massless tensor modes and the scalar mode, and considering that the GW is travelling along the z-axis with the speed of light $c=1$, we get the solution for $h_{\mu\nu}$ as $$\label{eq37}
 h_{\mu\nu}= \bar{h}_{\mu\nu}(t-z)-\bar{g}_{\mu\nu}\,\delta\phi(t-z).$$ Hence, for this solution we may take $q_{\mu} = \Omega\,(1,0,0,1)$ and $p_\mu = \omega\,(1,0,0,1)$. It should be noted that for $n\neq 2$ and $n\neq 1$ the above solution will contain a mixed state of massive scalar mode and breathing mode besides the tensor modes of polarization [@liang_polarizations_2017].

Geodesic Deviation and Polarization Modes
=========================================

To see the geodesic deviation, we substitute the solution (\[eq36\]) in the equation (\[eq34\]), which gives, $$\label{eq38}
R_{\mu\nu} = -\, \dfrac{1}{4}\, \bar{g}_{\mu\nu}\,\delta\phi - p_\mu p_\nu.$$ For the wave travelling along the $z$ direction, the nonzero components of $R_{\mu\nu}$ are $R_{tt}$, $R_{tz}$ and $R_{zz}$. Now to linear order we may write, $$R_{\mu\nu\alpha\beta} \approx \dfrac{1}{2} ( h_{\nu\alpha,\mu\beta}
+h_{\mu\beta,\nu\alpha} - h_{\mu\alpha,\nu\beta} - h_{\nu\beta,\mu\alpha}).$$ This can be further simplified as $$\label{eq39}
 R_{itjt} = -\, \dfrac{1}{2}(h_{ij,tt} + h_{tt, ij}).$$ Thus, using only the scalar part of the equation (\[eq37\]) we can write the geodesic equation as $$\label{eq40}
\ddot{x}^i = \dfrac{1}{2}\,\Big(\delta^i_j \delta\ddot{\phi} - \delta\phi_{\,,j}^{\,,i}\Big)x^j,$$ which gives, $$\begin{aligned}
\label{eq41}
 \ddot{x}& = \frac{1}{2}\,\delta \ddot{\phi} x,\\  
 \label{eq42}
 \ddot{y}& = \frac{1}{2}\,\delta \ddot{\phi} y,\\  
 \ddot{z}& = 0.\end{aligned}$$ Thus the geodesic equations show that there is no longitudinal component of the GW polarization for the pure $R^2$ model. So, for this model besides the tensor modes there exists only a scalar mode which is massless and pure in nature, known as the breathing mode. Again, from the equation (\[eq37\]) considering only the scalar part we may write, $$\label{eq43}
\delta \ddot{\phi}=-\,\omega^{2}\delta\phi.$$ Using this equation in the equation (\[eq41\]) we have, $$\label{eq44}
\ddot{x}+\frac{1}{2}\,\omega^{2}\delta\phi\, x = 0,$$ which for a GW propagating along $z$ axis takes the form: $$\label{eq45}
\ddot{x}+\frac{1}{2}\,\phi_{0}\,\omega^{2}\cos\omega(t-z) x = 0.$$ The solution of this equation (\[eq45\]) gives the time variation of deviation of $x$ at some fixed $z$. This equation is a special form of the well known Mathieu’s equation [@canosa_numerical_1971; @gheorghe_charged_0000] and the solution is graphically shown in Fig.\[fig2\] for two different set of parameters. It is seen that the time period of geodesic deviation depends highly on $\omega$. The results are identical along the $y$ direction also as clear from the equation (\[eq42\]).

![Geodesic deviation along $x$ for a fixed $z$. The top plot is for $\omega = 1.5$, $\phi_0  = 0.02$ and $z= 1$, and the bottom is for $\omega = 1$, $ \phi_0  = 0.2$ and $z = 1$. All parameters are in arbitary units.[]{data-label="fig2"}](01.eps)

![Geodesic deviation along $x$ for a fixed $z$. The top plot is for $\omega = 1.5$, $\phi_0  = 0.02$ and $z= 1$, and the bottom is for $\omega = 1$, $ \phi_0  = 0.2$ and $z = 1$. All parameters are in arbitary units.[]{data-label="fig2"}](02.eps)

Polarization Modes With Newman-Penrose Formalism
================================================

The Newman-Penrose (NP) formalism [@newman_approach_1962] can be used to find out the different polarization modes of GWs in a model. However, one major drawback of NP formalism is that it is only applicable to null waves. In $f(R)$ theory metric formalism, usually GWs have massive longitudinal mode of polarization due to which the NP formalism fails and shows deviated results [@liang_polarizations_2017]. But in the particular case of pure $R^2$ model, as we have already seen that the scalar field is massless and thus this allows us to use the NP formalism in the study of polarization modes of GWs in the theory. In NP formalism a GW is described with the help of six amplitudes $\lbrace \psi_2, \psi_3, \psi_4, \phi_{22} \rbrace$ representing six polarization modes in a particular coordinate system or frame [@newman_approach_1962]. All these amplitudes are defined as [@eardley_gravitational-wave_1973]: $$\begin{aligned}
\notag
\psi_2 &= -\,\dfrac{1}{6} R_{ztzt},\\\notag
\psi_3 &= -\,\dfrac{1}{2}R_{xtzt}+\dfrac{1}{2}i R_{ytzt},\\\notag
\psi_4 &= -\,R_{xtxt} + R_{ytyt}+ 2i R_{xtyt},\\\notag
\phi_{22} &= -\, R_{xtxt} - R_{ytyt}.\end{aligned}$$ Each of the complex amplitudes $\psi_3$ and $ \psi_4 $ are actually equivalent to two real amplitudes [@eardley_gravitational-wave_1973]. In Brans-Dicke theory, the massless scalar field appears as the breathing mode showing $ \phi_{22} = - \,R_{xtxt} -\, R_{ytyt}\neq 0.$ For our model these amplitudes are found as $$\begin{aligned}
 \notag
\psi_2& = 0,\\\notag
\psi_3 &= 0, \\\notag
\psi_4 &= \ddot{\bar{h}}_{xx}-\ddot{\bar{h}}_{yy}+ i\,(2\,\ddot{\bar{h}}_{xy}),\\\notag
\phi_{22} &= -\,2\,\delta\ddot{\phi}. \end{aligned}$$ These results show that our model, i.e. pure $R^2$ model gives three modes of massless polarization: two for tensor modes and one for the breathing mode as found earlier, and according to Lorentz - invariant $E(2)$ classification of plane waves, the model results GWs of class $N_3$ and all modes are independent of the observer [@eardley_gravitational-wave_1973].

Pulsar Timing Arrays and Correlation of Polarization Modes
==========================================================

Pulsar timing arrays (PTAs) can be used to detect the polarization modes of GWs. A good number of works have been going on in this field and PTAs are found to be effective in the detection of extra polarization modes to be present in GWs. Hence they can be used as a tool to test modified theories of gravity [@lee_pulsar_2008; @mingarelli_characterizing_2013; @lee_detecting_2010].

Presence of GWs disturb the null geodesic of the signals from pulsars. Due to this reason the time of arrival of GWs changes. So, by tracing the changes in time of arrival of the radio signals from the pulsars, it is possible to detect the GWs. For the mathematical treatment of the PTA procedure let us consider a PTA detecting radio signals in the regime of GWs. In case of GWs from the pure $R^2$ model, the information about the source is carried by the three polarization amplitudes, viz., $h_+(t)$, $h_\times(t)$ and $h_\phi(t)$, where the first two stands for the Einstein or tensor modes and the third one for scalar (breathing) mode of polarization. Thus the GW signal from a source for the pure $R^2$ model can be wriiten as $$\label{eq46}
h_{ij}(t,\hat{\Omega}) = e^+_{ij}(\hat{\Omega})h_+(t,\hat{\Omega}) + e^\times_{ij}(\hat{\Omega})h_\times(t,\hat{\Omega}) + h_\phi(t,\hat{\Omega}).$$ Here $e^+_{ij}$ and $e^\times_{ij}$ are the polarization tensors as given by, $$\begin{aligned}
\nonumber
e^+_{ij}(\hat{\Omega}) = \hat{m}_i \hat{m}_j - \hat{n}_i \hat{n}_j,\\\nonumber
e^\times_{ij}(\hat{\Omega}) = \hat{m}_i \hat{n}_j + \hat{n}_i\hat{m}_j.\end{aligned}$$ Hence, these polarization tensors are uniquely defined once the unit vectors $\hat{m}$ and $\hat{n}$ used to describe the principal axies of wave are specified. $\hat{\Omega}$ is the direction of GW propagation given by $\hat{m} \times \hat{n}$. As per our previous convention, we set $\hat{\Omega} = \hat{z}$, and consider a null vector $S^\mu$ that points to the Solar system barycentre from the pulsar in Minkowski spacetime. In purturbed spacetime vector $S^\mu$ will change to another vector $\sigma^\mu$ as given by, $$\label{eq47} 
\sigma^\mu = S^\mu - \dfrac{1}{2}\eta^{\mu\nu}h_{\mu\nu}S^\nu + \delta \phi\, \eta^{\mu\nu}b_{\mu\nu} S^\nu,$$ where the second part of this equation is due to tensor modes of purturbation and the third part is due to the scalar mode of purturnation in spacetime with $b_{\mu\nu}$ is a unit breathing mode matrix having two non-zero unit components $b_{11}$ and $b_{22}$, obtained by the application of transverse condition to this mode of GWs. It is to be noted that for this mode we can not apply the traceless condition as the application of this condition to this mode will retain only the tensor modes or GR modes [@kausar_gravitational_2016]. Now, If we define $S^\mu$ as $S^\mu = \nu(1,-\alpha,-\beta,-\gamma)$, where $\alpha,\beta$ and $\gamma$ are direction cosines, and $\nu$ is the frequencey of the radio pulses from the source, then from equation (\[eq47\]), we may write, $$\begin{aligned}
\notag 
\sigma^t& = \nu\,\\\notag
\sigma^x &= -\,\nu\,\big[\alpha\,\big(1-\dfrac{1}{2}h_+ - \delta\phi\big)-\dfrac{\beta}{2}h_\times\big],\\\notag
\sigma^y& = -\,\nu\,\big[\beta\big(1+\dfrac{1}{2}h_+ + \delta\phi\big)-\dfrac{\alpha}{2}h_\times\big],\\\notag
\mbox{and}\;\;
\sigma^z &= -\,\nu\,\gamma.\end{aligned}$$ The radio pulses from the pulsar follow a null geodesic through spacetime. The geodesic equation of the pulses with the affine parameter $\lambda$ is $$\label{eq48}
\dfrac{d\nu}{d\lambda} = -\,\nu^2\alpha\,\beta\,\dot{h}_\times + \dfrac{1}{2}\,\nu^2 \big[\beta^2(\dot{h}_+ - \delta\dot{\phi}) - 2 \alpha^2 (\dot{h}_+ + \delta\dot{\phi}) \big].$$ Using this equation in the geodesic deviation equation we can have the frequency shift at Solar system barycentre as given by, $$\begin{aligned}
\label{eq49}
z(t, \hat{\Omega}) &=\; \mid \dfrac{\nu(t)- \nu_0}{\nu_0} \mid \\\notag
&= \dfrac{1}{2(1+\gamma)}\Big[\alpha^2 (\Delta h_+ + \Delta \delta\phi)
 - \beta^2(\Delta h_+ - \Delta \delta \phi)\Big] +\dfrac{\alpha \beta}{(1+\gamma)}
 \Delta h_\times,\end{aligned}$$ where $\nu(t)$ is the frequency observed at Solar system barycentre. Thus it is seen that the presence of GWs will give rise to a frequency shift and these shifts can be observed with the help of PTAs. It is also clear that the presence of massless breathing mode results in a contribution of a shift in the frequency besides a contribution from the usual tensor (GR) modes. Following [@lee_pulsar_2008], for stochastic GW background with normalized frequency $f$, the correlation function for GR modes and breathing mode of polarization are calculated. This function for tensor modes is found as $$\label{eq50}
C^{+,\times}(\theta) = \xi^{GR}(\theta)\int_0^\infty\dfrac{|h_c^{+,\times}|^2}{24\pi^2 f^3}\,df,$$ where $$\xi^{GR}(\theta)=\dfrac{3\,(1-\cos\theta)}{4}\log\,(\dfrac{1-\cos\theta}{2}) + \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 $$\label{eq51}
C^{b}(\theta) = \xi^{b}(\theta)\int_0^\infty\dfrac{|h_c^{b}|^2}{12\pi^2 f^3}\,df,$$ where $$\xi^{b}(\theta) = \dfrac{1}{8}\big[\cos\theta + 3 + 4\,\delta(\theta)\big].$$ The correlation coefficients $\xi^{GR}$ and $\xi^{b}$ as a function of $\theta$ are plotted in the Fig.\[fig3\] It is seen that provided the cases in [@lee_pulsar_2008] hold good, PTAs can effectively distinguish between the breathing mode and tensor modes present in GWs. It needs to mention that the correlation versus angular separation curve for the tensor modes of polarizations was first obtained by Hellings and Downs in 1983 and hence it is usually known as Hellings-Downs (HD) curve [@hellings_upper_1983; @jenet_understanding_2014].

![Variation of correlation coefficient $\xi(\theta)$ of different polarization modes with respect to angular separation between two pulsars.[]{data-label="fig3"}](correlation_d.eps)

Summary and Conclusions
=======================

In this work, we have used the $f(R)$ gravity power law model to study the polarization modes of Gravitational Waves in de Sitter spacetime. We have seen that the field equation in de Sitter spacetime gives a massive scalar mode with a mass term $m^2 = \dfrac{(2-n)\bar{R}}{(n-1)}.$ The mass term varies widely with the exponent term $n$ and the background curvature or the de Sitter curvature, and the mass of the scalar field becomes zero when background curvature is zero or $n=2$. Later using the stability condition for the theory in de Sitter spacetime we have seen that for a constant curvature $\bar{R}$ the theory has stable solutions for $n=2$. It has been observed that for this particular case of the model, the massive longitudinal mode of polarization vanishes. Thus, this is the only case in $f(R)$ theory in metric formalism treatment where only three polarization modes are obtained and the third mode is a pure breathing mode. To validate this result we’ve studied the geodesic deviations for the pure $R^2$ model explicitly and later we used NP formalism to confirm the validity of our result. The absence of the massive scalar field in this model allows us to use NP formalism, which is a powerful tool to check the polarization contents of null GWs. The results from the analysis show absence of massive polarization mode of the GWs in the theory, which establish that the polarization modes and the mass of the scalar field in $f(R)$ gravity are model dependent. Another important result obtained from this work is that the scalar field is independent of the background curvature in pure $R^2$ model of $f(R)$ gravity and evidently chameleonic behaviour is not observed here. The absence of massive longitudinal mode makes this pure $R^2$ model different from the other $f(R)$ gravity models and hence a more detail study in this model is required.

The pure $R^2$ model can be extended by adding other terms to the action which can account for the missing part giving rise to mixed polarization states of GWs besides tensor modes of polarization [@kausar_gravitational_2016; @liang_polarizations_2017]. Those extended models can be constrained with the experimental results obtained so far. It is expected that future experiments can provide better constraints to the $f(R)$ gravity models using which the existence of extra polarization modes can be hopefully confirmed, which in turn give options to test the reliability of $f(R)$ gravity in the modifications and extensions of GR. Extensions like nonminimal matter field coupling and other modified gravity theories can be included as the future aspects of this type of works, which can be tested with future experiments on GWs.

Spin projection operators
=========================

The explicit form of the spin projection operators that span the solutions space to the linearized field equations are given below [@van_nieuwenhuizen_ghost-free_1973; @van_nieuwenhuizen_supergravity_1981]: $$\begin{aligned}
\notag
P^2&=\frac12(\theta_{\mu\rho}\theta_{\nu\sigma}+\theta_{\mu\sigma}\theta_{\nu\rho})-\frac13\theta_{\mu\nu}\theta_{\rho\sigma},\\\notag
P^1_e&=\frac12(\theta_{\mu\rho}\omega_{\nu\sigma}\theta_{\mu\sigma}\omega_{\nu\rho}-\theta_{\nu\rho}\omega_{\mu\sigma}+\theta_{\nu\sigma}\omega_{\mu\rho}),\end{aligned}$$ $$\begin{aligned}
\notag
P^1_b&=\frac12(\theta_{\mu\rho}\theta_{\nu\sigma}-\theta_{\mu\sigma}\theta_{\nu\rho}),\\\notag
P^1_m&=\frac12(\theta_{\mu\rho}\omega_{\nu\sigma}+\theta_{\mu\sigma}\omega_{\nu\rho}+\theta_{\nu\rho}\omega_{\mu\sigma}+\theta_{\nu\sigma}\omega_{\mu\rho}),\\\notag
P^1_{em}&=\frac12(\theta_{\mu\rho}\omega_{\nu\sigma}+\theta_{\mu\sigma}\omega_{\nu\rho}-\theta_{\nu\rho}\omega_{\mu\sigma}-\theta_{\nu\sigma}\omega_{\mu\rho}),\\\notag
P^1_{me}&=\frac12(\theta_{\mu\rho}\omega_{\nu\sigma}-\theta_{\mu\sigma}\omega_{\nu\rho}+\theta_{\nu\rho}\omega_{\mu\sigma}-\theta_{\nu\sigma}\omega_{\mu\rho}),\\\notag
P^0_{sw}&=\frac1{\sqrt3}\theta_{\mu\nu}\omega_{\rho\sigma},\\\notag
P^0_{ws}&=\frac1{\sqrt3}\omega_{\mu\nu}\theta_{\rho\sigma},\\\notag
P^0_s&=\frac13\theta_{\mu\nu}\theta_{\rho\sigma},\\\notag
P^0_w&=\omega_{\mu\nu}\omega_{\rho\sigma},\end{aligned}$$ where the transverse and longitudinal projection operators in the momentum space are respectively $$\theta_{\mu\nu}=\delta_{\mu\nu}-\frac{\nabla_\mu \nabla_\nu}{\square},\qquad \omega_{\mu\nu}=\frac{\nabla_\mu \nabla_\nu}{\square}.$$

[99]{} H. Weyl, Annalen der Physik [**364**]{}, 101 (1919). A. S. Eddington, *The mathematical theory of relativity*, Cambridge Univ. Press, Cambridge (1924). R. Utiyama and B. S. DeWitt, J. Math. Phys. [**3**]{}, 608 (1962).

T.P. Sotiriou and V. Faraoni, Rev. Mod. Phys. [**82**]{}, 451 (2010) \[arXiv:0805.1726\]. B. P. Abbott and R. Abbott et al., Phys. Rev. Lett. [**116**]{}, 061102 (2016) \[doi: 10.1103/PhysRevLett.116.061102\]. B. P. Abbott and R. Abbott et al., Phys. Rev. Lett. [**116**]{}, 241103 (2016) \[doi: 10.1103/PhysRevLett.116.241103\]. B. P. Abbott and R. Abbott et al., Phys. Rev. Lett. [**118**]{}, 221101 (2017) \[doi: 10.1103/PhysRevLett.118.221101\]. B. P. Abbott and R. Abbott et al., Phys. Rev. Lett. [**119**]{}, 141101 (2017) \[doi: 10.1103/PhysRevLett.119.141101\]. B. P. Abbott and R. Abbott et al., Phys. Rev. Lett. [**119**]{}, 161101 (2017) \[doi:10.1103/PhysRevLett.119.161101\]. Y. S. Myung, arXiv:1608.01764 (2016). D. Liang, Y. Gong, S. Hou, and Y. Liu, Phys. Rev. D [**95**]{}, 104034 (2017) \[arXiv:1701.05998\]. L. G. Jaime, L. Patio, and M. Salgado, Phys. Rev. D [**87**]{}, 024029 (2013) \[arXiv:1212.2604\]. L. Amendola, D. Polarski, and S. Tsujikawa, IJMP D [**16**]{}, 1555 (2007). L. Amendola, D. Polarski, and S. Tsujikawa, Phys. Rev. Lett. [**98**]{} 131302 (2007) \[arXiv:astro-ph/0603703\]. L. Amendola, R. Gannouji, D. Polarski, and S. Tsujikawa, Phys. Rev. D [**75**]{}, 083504 (2007) \[arXiv:gr-qc/0612180\]. S. Capozziello, S. Carloni, and A. Troisi, arXiv:astro-ph/0303041 (2003). S. Carloni, P. K. S. Dunsby, S. Capozziello, and A. Troisi, Class. Quant. Grav. [**22**]{}, 4839 (2005) \[arXiv:gr-qc/0410046\]. U. D. Goswami and K. Deka, IJMP D [**22**]{}, 1350083 (2013) \[arXiv:1303.5868\]. V. Faraoni and S. Nadeau, Physical Review D [**72**]{}, 124005 (2005) \[arXiv:gr-qc/0511094\]. J. D. Barrow and A. C. Ottewill, J. Phys. A: Math. Gen. [**16**]{}, 2757 (1983).

V. Faraoni and S. Nadeau, Phys. Rev. D [**72**]{}, 124005 (2005) \[arXiv:gr-qc/0511094\].

A. Nunez and S. Solganik, arXiv:hep-th/0403159 (2004). P. van Nieuwenhuizen, Nucl. Phys. B [**60**]{}, 478 (1973). P. van Nieuwenhuizen, Phys. Rep. [**68**]{}, 189 (1981). T. Chiba, Phys. Lett. B [**575**]{}, 1 (2003). J. Canosa, J. Comput. Phys. [**7**]{}, 255 (1971). V. N. Gheorghe and G. Werth, *Charged particle traps: Physics and techniques of charged particle field confinement* (Physica-Verlag NetLibrary, Inc. distributor, Heidelberg; Boulder, 2005), ISBN 9783540265764.

E. Newman and R. Penrose, J. Math. Phys. [**3**]{}, 566 (1962). D. M. Eardley, D. L. Lee, and A. P. Lightman, Phys. Rev. D [**8**]{}, 3308 (1973). K. J. Lee, F. A. Jenet, and R. H. Price, The Astrophysical Journal [**685**]{}, 1304 (2008). C. M. F. Mingarelli, T. Sidery, I. Mandel, and A. Vecchio, Phys. Rev. D [**88**]{}, 062005 (2013) \[arXiv:1306.5394\]. K. Lee, F. A. Jenet, R. H. Price, N. Wex, and M. Kramer, The Astrophysical Journal [**722**]{}, 1589 (2010) \[arXiv:1008.2561\]. R. W. Hellings and G. S. Downs, The Astrophysical Journal [**265**]{}, L39 (1983). F. A. Jenet and J. D. Romano, arXiv:1412.1142 (2014). H. R. Kausar, L. Philippoz, and P. Jetzer, Phys. Rev. D [**93**]{}, 124071 (2016) \[arXiv:1606.07000\].
