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

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 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 are two extra scalar modes of polarization of GWs present in the theory. These two extra polarization modes present in a mixed state, out 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 massive longitudinal 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.


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 [1][2][3][4]. 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 [5]. 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 [6] 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 [7], Hu-Sawicki model [8], Tsujikawa model [9][10][11][12] 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 [5,13]. It is found that total number of polarization modes of GWs that can exist in f (R) theories of gravity is 4 [14].
Recent studies show that there exist a massless breathing mode and a massive longitudinal mode of polarization of GWs in f (R) theories of gravity. These two extra polarization modes of GWs exist in a mixed state and their existence can be checked with Newmann-Penrose formalism and also by geodesic deviation equation in the theory [15]. However, another study shows that these two extra polarization modes are model dependent and there exists a model in f (R) gravity where the massive longitudinal mode of polarization vanishes [16].
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, the model has been constrained. In section four, we have studied the polarization modes of GWs present in the model by using the perturbed field equation and the Newmann-Penrose formalism. In the fifth 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. 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. Thus a new f (R) gravity model capable of explaining the drawbacks of GR might have a complete new set of physics 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 free 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 free parameters, over and above to the existing ones, as given by: where α and β 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: The first correction factor is estimated by two free parameters α and R c . Similarly, the second correction factor has also two free parameters β and R c and it mimicks 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 [17][18][19][20]: (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 de-Sitter solution for a model can be obtained by solving the equation, where R 0 is the scalar curvature of present time. For our model this equation takes the form: where x = R 0 /R c and R c = 0. The value of R 0 satisfying the above complex equation is the stable de-Sitter background curvature. If we pick a particular case say R c = R 0 , then the de-Sitter solution is R 0 = 0, i.e. Minkowski spacetime.
(3) A positive effective gravitational coupling, leading to f ′ (R) > 0. By putting our model in this condition gives, (4) A stable cosmological perturbation and a positivity of the GWs for the scalar mode, causing to f ′′ (R) > 0. Using our model in this condition we find, This condition ensures the absence of tachyonic instabilities in the model, i.e. this condition results m 2 φ > 0. (5) An asymptotic behaviour to the ΛCDM model in the large curvature region. In case of our model, at large curvature region i.e. R → ∞, we have f (R) − R → − α/2 − β = constant, which mimicks the ΛCDM model in the large curvature region. Again, at 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.

A. Scalar Degrees of Freedom in Jordan Frame
The action of a generic f (R) gravity model is given as In the above equation, the function f (R) stands for any arbitrary function of Ricci curvature scalar R, g µν is the metric, κ 2 = 8πG = M −2 pl and = c = 1. Here M pl ≈ 2 × 10 18 GeV is the reduced Planck's mass. L m g µν ,ψ is the Lagrangian for a matter fieldψ. Variation of the action (2) with respect to the metric gives the following field equation: where δg µν is the matter energy-momentum tensor and f ′ (R) = ∂ R f (R).

Trace of Eq. (3) is
where T = g µν T µν is the trace of T µν . 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 then for our model the field φ becomes, This shows that the scalar curvature R can be expressed as a function of the field φ. From this definition of scalar field φ, we may view the trace Eq. (4) as an effective scalar field equation of Klein-Gordon with the following identification [21]: where V is the potential of the scalar field φ. Thus, the trace Eq. (4) can be written as a Klein-Gordon type equation for the scalar field φ as given by, where V ef f is effective potential of the field and is define as At far away from the source or in absence of any matter source V ef f ≡ V . Again, from the stationary condition: we can have φ = φ 0 satisfying φ 0 = f ′ (R 0 ). From this condition, the mass of the scalar field (or the scalaron mass) can be obtained by differentiating the Eq. (9) with respect to φ as Here R 0 is the background curvature corresponding to φ 0 . From the above equation, we can see that avoiding of tachyonic instabilities demands f ′ (R0) f ′′ (R0) − R 0 ≥ 0 and to keep the mass term finite we need f ′′ (R 0 ) = 0. For our model the mass term is found as For R 0 = 0, this equation gives 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 φ of the scalar field depends on the model parameters α, β and R c . Fig. 1 shows the variation of m 2 φ 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 α takes value closer to parameter β. By increasing the difference between β and α (i.e. for β ≫ α) the hump can be minimized. An increase of α 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.

Scalar Tensor Equivalence of the Model
To see the origin of the scalar field in the theory, we would like to rewrite the action (2) by introducing a new auxiliary scalar field ψ as [15,22] where f ψ = ∂f (ψ) ∂ψ . Now, varying this equation with respect to the new auxiliary scalar field ψ we get, For finiteness of the previously defined mass square term of the scalar field (see Eq. (11)), we have f ′′ (R) = 0, which is in terms of ψ, f ψψ (ψ) = 0. With this condition the above equation gives, R = ψ. Substituting of this result in the action (14) we can recover the original action (2). Moreover, the quantum stability condition demands that f ′′ (R) ≥ 0. This along with the finiteness condition of the mass of the scalar field demands that f ′′ (R) > 0. Thus, it is always possible to have a scalar tensor representation of f (R) theory of gravity. Redefining the previously defined scalar field φ in terms of the new auxiliary field ψ as the action (14) can be rewritten as where is the potential of the scalar field. To be precise, this term f ψ (ψ) ψ − f (ψ) 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, For f (R) = R and hence for f ′ (R) = 1, the action (17) recovers GR giving the potential term V = 0. Using Eq.
(1) in this Eq. (18), the scalar field potential for our model can be obtained as The variation of the potential V with respect to R 0 for different values of α and β is shown in Fig 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 α and β. 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 α moves closer to the value of β. In fact, this dip in the potential is responsible for the hump in the mass square curve for the corresponding values of α and β. This dip region of the potential curve almost eliminates in the case when α ≪ β. Moreover, with the increasing values of both α and β, the potential comparatively increases after the dip region or without the dip region.

B. 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 [23,24]:g µν = f ′ (R)g µν , which for our model takes the form:g Consequently, in this frame with L m = 0, the action changes to [23,24] where the scalar field and V (φ E ) is potential of the scalar field in this frame, and is given by, The Eq. (22) shows the dependency of the Einstein frame scalar field φ E on the scalar curvature R. Using our model, this potential (23) can be expressed as where x = R 0 /R c and χ = x 4 + 1. Hence, the mass square term of the scalar field in the Einstein frame is In Minkowski space i.e. for R 0 = 0, this equation takes the form: From Eq.s (24) and (25), 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 φ E as a function of background curvature is shown. The field φ E at R 0 = 0 and R 0 → ∞ is independent of α, whereas it is independent of β only at R 0 → ∞. Moreover, φ E is non-zero at R 0 = 0 and tends to zero at R 0 → ∞, which is obvious from it's expression. Again from Eq.s (6) and (22) it is clear that all these behaviours of φ E should be applicable to φ also, but with positive values of φ for all values of R 0 . Variation of the potential (24) as a function of φ E is shown in Fig. 6. It is seen from this figure that the minimum of the potential moves towards the higher value of φ E when α ≪ β than the case when α ∼ β. Obviously, similar behaviour can be attributed to the potential (19) as a function of the field φ. 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.

C. 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 φ = f ′ (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 [25]. According to Guo, a model can pass solar system tests if it satisfies the following conditions: We've calculated the above functions numerically for our model (see Table I). 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 free parameter R c . However, a simple and effective way to make the model solar system viable is to decrease all the free parameters (i.e. α, β and R c ) sufficiently (see Table I). 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.  to constrain an f (R) gravity model [26][27][28][29][30][31][32]. 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 [26,27] and [29]. In Ref. [26], 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 σ 8 and the matter mass-energy density Ω m from Planck Sunyaev-Zeldovich (PSZ) cluster number counts and also from the CFHTLens weak lensing tomography measurements. Combining CMB, BAO and σ 8 − Ω m relationship from the PSZ catalog [33], they obtained a bound which is still better than the bounds obtained from the GW event GW170817 [29]. The bound on the parameter f ′ (R) reported by them is with 95% confidence level at upper bound [26]. On the other hand in Ref. [27] a constraint was introduced on the Compton wavelength λ g of the graviton. From their study, we have a constraint on λ −1 g as given by, with 90% confidence level on upper bound [27]. Now, we have computed the values of f ′ (R) − 1 and λ −1 g for our model by taking into consideration of above cited respective upper bounds with the corresponding confidence levels to constraint our free model parameters R c , α and β. The results of this computation along with the contour plots are shown in Fig. 8. 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 ′ (R) − 1 and 90% confidence level for λ −1 g (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 ′ (R) − 1 and λ −1 g and any value lower than the boundary value is viable. This point corresponds to the galaxy cluster Abell 1689 data [27]. In the plots we have considered three sets of the free parameters and we see that the smaller values of α allow the model to pass the constraints easily. In the Fig. 8, we have also shown 3 other points corresponding to galaxy clusters Abell 262, Abell 1991 and Abell 383 data from the Ref. [28]. All these points lie within the confidence level contours.
The model can be constrained by using the GWs event GW170817 also. In a recent study [29], f (R) gravity was constrained by using the GW170817. They provided a bound on f ′ (R), which is Now for the sack of simplicity, we would like to set the arbitrary characteristic constant R c = R 0 , where as earlier R 0 is the background curvature. Under this assumption the bound (32) for our model takes the following form: .
This is the constraint on the model parameters provided by the GWs event GW170817. Assuming α and β both to be positive quantities, we can simplify the above constraint to the following form: It should be noted that, here the constraint range of α is β dependent, which is constraint by the relation on the left side of the logical AND operator ∧.  [27] (the larger contour) and with 68% confidence level on the both (the smaller contour). The central red dot denotes the λg corresponding to galaxy cluster Abell 1689 [27], blue dot corresponds to Abell 262, black one corresponds to Abell 1991 and yellow one corresponds to Abell 383 [28]. data  [26] and λg [27] .
In Fig. 9, we have shown the variations of λ g with respect to R 0 , f ′ (R)−1 with respect to R 0 and m g with respect to f ′ (R)−1 for the values of the model parameters used in the contour plots in Fig. 8. The hump in the mass of the scalar field encountered in the Fig. 1 and in the Fig. 6 are also present in the λ g vs. R 0 curves. However, as mentioned earlier, this hump vanishes when (β − α) ≫ 0. For higher curvatures, λ g rapidly moves towards zero. The function f ′ (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 ′ (R) − 1 mass of the scalar field decreases initially at a faster rate and then becomes almost constant at later stage.

IV. 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.

A. 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ḡ µν we may express the spacetime metric to the first order of perturbation value h µν , which is usually usually very small, as Now, expanding the Ricci tensor and the Ricci scalar upto the first order of h µν , we may write: and Similarly, we may write for the f (R) and f ′ (R) as whereR is some constant curvature. Thus, due to the perturbation in spacetime the trace equation (4) can be rewritten as where we have used T µν = 0 for the empty space or far away from the source. Fixing the gauge to be harmonic gauge with ∇ µ h µ ν = 1 2 ∇ ν h, which after operating by ∇ ν we find, An important point to be mentioned here is that, the Eq. (3) is also satisfied by another solution: R µν = Λg µν =R µν , giving This equation actually corresponds to the Eq. (10), which is the stationary condition used earlier in the Sec. II. In empty space, this equation has the form: and it leads to have the equation, This is the stationary condition (10) of the scalar field potential in the empty space corresponding to the constant curvatureR of spacetime. Using the Eq.s (37), (41) and (44) in Eq. (40), we get Now, we would like to define h = m 2 h, where m is the mass of the associated scalar field. Using this definition in the above equation we obtain, This is a quadratic equation in m 2 and solution for m 2 gives, and 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. (11). Thus the term m 2 given by Eq. (47) is exactly same as the scalar field mass square term m 2 φ given in Eq. (12) for our model, whenR = 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ḡ µν = η µν , i.e. the Minkowski metric and the background curvatureR = 0. In this case, the Ricci scalar slowly varies near zero, i.e. R ≃ 0 + δR. Hence, for the Minkowski space the Eq. (35) can be written as And to the first order of h µν , we get where h = η µν h µν . So for our model, to the first order of perturbation, the Eq. (3) becomes Taking the trace of this Eq. (52), we get where with α > 0 and β > 0. This is also exactly the same mass square term m 2 φ R0 = 0 in Minkowski space given by the Eq. (13) for our model. Next, we introduce a variableh The trace of this variable ish Using this Eq. (55) in the variable (54) we find, From Eq. (54) and Eq. (56), one can easily see that both h µν andh µν are interchangeable, i.e. replacingh µν by h µν and vice-versa in Eq. (54) gives Eq. (56). Again, under an infinitesimal coordinate transformation, x µ → x µ′ = x µ + ς µ , we have The trace of this equation is The trace trace of this equation gives,h Here we raise or lower the indices with the help of η µν , i.e. with the Minkowski metric. The Lorentz gauge condition ∂ µh′ µν = 0 can be obtained if ς µ satisfies ς ν = ∂ µh µν . 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. ∂ µh µν = 0 andh = η µνh µν = 0 [15,[34][35][36]. By using the transverse traceless gauge condition and substituting the Eq. (56) into Eq. (50), we get Plugging Eq. (61) into Eq. (52), we obtain Combining Eq.s (53) and (62), we get which is the wave equation of the massless tensor field. The solution to this Eq. (63) is [34,35] h µν = e µν exp(iq µ x µ ) + c.c., where η µν q µ q ν = 0 and q µ e µν = 0. Whereas the solution to the massive scalar field Eq. (53) is given by [34,35], where η µν p µ p ν = − m 2 0 . Assuming the GWs propagation direction along z, the general solution can be written as [15] h µν =h µν (t − z) whereh µν is transverse and traceless and it represents the standard spin-2 graviton, and ψ 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. ψ(vt − z) can be expressed as and hence the mass of the field in terms of k and ω is

B. Calculation of Newmann-Penrose Quantities of the Model
In 1973, a powerful method was introduced in Ref. [37] 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 [38]. In the Newmann-Penrose formalism, there are ten Ψ's, nine Φ's, and a Λ, which are algebraically independent and represent the irreducible parts of the Riemann tensor R λµκν . They are known as Newmann-Penrose quantities. But in case of plane waves or nearly plane waves, the differential and symmetry properties of R λµκν reduce the number of independent, nonvanishing components, to six. Hence, in this formalism, the set {Ψ 2 , Ψ 3 , Ψ 4 , Φ 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 [37]: It should be noted that, Ψ 3 and Ψ 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 polarization amplitudes for our model as where m 2 0 is given by Eq. (68) and C 1 = 2α−πβ π(1−β)Rc (see Eq. (66)). From the above expressions we can calculate the Newmann-Penrose quantities. Note that, above results suggest, there are 4 non-zero polarization amplitudes in the theory. Using Eq.s (80), we've calculated the Newmann-Penrose quantities for the model as Thus the E(2) classification of the model is II 6 . The model exhibits 4 polarization modes viz., tensor plus, tensor cross, scalar transverse massless breathing mode and scalar longitudinal massive mode of polarization. The breathing mode and the longitudinal mode exist in a mixed state. If m 0 = 0, the massive longitudinal mode will vanish, giving Ψ 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 Ψ 2 = 0 and Φ 22 = 0. When both m 0 and ψ(vt − z) vanish, the theory reduces to GR giving only tensor modes of polarizations.

V. 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 [40]. Later in 1982, the first millisecond pulsar was discovered [41]. 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 [42,43]. After three years, the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) in North America was founded [44,45]. 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 [46]. 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 [47,48]. In this study, we have used some selected data from PPTA [42] and IPTA [47,48].
From the Ref. [49], 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 [49] C +,× (θ) = ξ GR (θ) where and θ is the angular separation between two pulsars. For the scalar modes it is [49] where ξ b (θ) = 1 8 cos θ + 3 + 4 δ(θ) .
The normalized correlation function in general is given by: 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. [50][51][52] in which the timing residual induced by GWs is expressed as here S = 2 (1 + (c/ω g )k g ·n) gives the dispersion relation of the GWs, the terms ω g and k g connect the mass of the longitudinal mode of polarization m g by the relation m 2 g = ω 2 g − k 2 g , A ij ≡n inj and H ij = τ 0 h ij (τ, 0) − h ij (τ − |D|/c, D) dτ . D is the displacement vector from the observer to the pulsar andn i andn j are two unit vectors pointing to two pulsars. The correlation coefficient C between two different pulsars is given by where the sub-scripts are indices for the pulsars. With these assumptions and following Ref. [49], the correlation functions are calculated numerically for different values of m g (see Fig. 10). In terms of our model, All these numerically calculated correlation functions are plotted with respect to θ in Fig. 10. In this figure we have used three values of m g for the following sets of parameters: 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 [42,47,48,53,54] as mentioned above, which are also plotted in Fig. 10. 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.

VI. 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 4 polarization modes of GWs viz., tensor plus mode, tensor cross mode, scalar breathing mode and scalar 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. These two non tensor and extra polarization modes exist in a mixed state along the scalar degrees of freedom. 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 σ 8 − Ω m relationship from the PSZ catalog and Abell 1689 galaxy cluster data. The model can be easily constrained with the help of the free parameters R c , α and β. Further, we have also constrained the model using the GW event GW170817 under the assumption R c = R 0 . 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 [14]. 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.