{'fraction_non_alphanumeric': 0.03557111686202872, 'fraction_numerical': 0.008651748109808698, 'mean_word_length': 4.718949506674405, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmology of F(R) nonlinear massive gravity

The theory of nonlinear massive gravity can be extended into the F(R) form as developed in Phys.Rev.D90, 064051 (2014). Being free of the Boulware-Deser ghost, such a construction has the additional advantage of exhibiting no linear instabilities around a cosmological background. We investigate various cosmological evolutions of a universe governed by this generalized massive gravitational theory. Specifically, under the Starobinsky ansantz, this model provides a unified description of the cosmological history, from early-time inflation to late-time self-acceleration. Moreover, under viable F(R) forms, the scenario leads to a very interesting dark-energy phenomenology, including the realization of the quintom scenario without any pathology. Finally, we provide a detailed analysis of the cosmological perturbations at linear order, as well as the Hamiltonian constraint analysis, in order to examine the physical degrees of freedom.


I. INTRODUCTION
Einstein's General Relativity (GR) is commonly acknowledged as the standard theory of gravitation at distance scales which are sufficiently large compared to the Planck length. As a foundation of modern cosmology, GR has been greatly successful in explaining the dynamics of the universe throughout the whole thermal expanding history. Furthermore, recent high-precision cosmological observations [1] strongly support the standard inflationary Big-Bang paradigm, in which our universe experienced a period of inflation at early times, then evolved through radiation and matter domination, and eventually entered a second cosmic acceleration today. Although such a picture can be described in the GR framework, with some additional scalar degrees of freedom, modifications of GR provide an alternative explanation.
Amongst all modified gravity theories, the F (R) paradigm is one simple theory that can sufficiently describe the properties of higher-order gravitational effects, by extending the gravitational Lagrangian as an arbitrary function of the Ricci scalar (see [2][3][4] for reviews). Hence, this scenario is expected to have important effects at high-energy scales, that were realized in the very early universe. In particular, the Starobinsky model [5] is one of the earliest inflationary models, and it lies in the center of the best-fit regime of the recently released observational results by the Planck collaboration [6]. Moreover, the F (R) modifications of GR can also be applied to explain the cosmic acceleration of the late-time universe [7], or inflationary and late-time epochs in a unified way [8].
On the other hand, the question whether the graviton can acquire a mass, has attracted the interest of theorists for decades. Initiated by Fierz and Pauli [9], the * Electronic address: yifucai@physics.mcgill.ca † Electronic address: Emmanuel˙Saridakis@baylor.edu subsequent full nonlinear formulation of massive gravity was found to suffer from a severe problem, namely the existence of the so-called Boulware-Deser (BD) ghost [10]. This fundamental problem puzzled physicists until recently, where a specific nonlinear extension of massive gravity was proposed by de Rham, Gabadadze and Tolley (dRGT) [11]. Within this new theory, through a Hamiltonian constraint analysis [12] and an effective field theory approach [13], it is proven that the BD ghost can be eliminated by a secondary Hamiltonian constraint (see [14] for a review). Apart from the theoretical interest, this nonlinear construction has the additional advantage that it can provide an explanation of the observational evidence of the late-time cosmic acceleration. In particular, the graviton potential can effectively mimic a cosmological constant by fine tuning a sufficiently small value for the graviton mass [15,16]. However, the cosmological perturbations around these solutions exhibit in general severe instabilities [17].
Having the above in mind, it is interesting to consider the possibility of combining the F (R) paradigm and the recent nonlinear massive gravity, as proposed in [18]. Such an extension allows for a significantly enhanced class of cosmological evolutions, especially describing inflation and late-time acceleration in a unified way. However, the important advantage is that this construction exhibits stable cosmological perturbations at the linear level. In the present work we perform a detailed analysis on the model of F (R) massive gravity as an accompanied work to [18]. We first perform a detailed generalized Hamiltonian constraint analysis in order to show the absence of the BD ghost. Moreover, we investigate in detail various cosmological evolutions, from inflation to dark energy eras, using viable F (R) ansantzes and focusing on observables such as the dark-energy equation-of-state parameter. Finally, we perform a detailed analysis of the linear cosmological perturbations.
The paper is organized as follows. In Section II we briefly review the paradigm of F (R) nonlinear massive gravity. We then search for cosmological solutions in both Einstein and Jordan frames in Section III. In Section IV we present some examples of the interplay between the F (R) and the massive sectors of the gravitational action, since this construction allows for a large variety of cosmological evolutions in agreement with the observed behavior. In Section V we perform a perturbation analysis at linear order around the cosmological background, which reveals only a single propagating scalar mode. Since the kinetic term of the scalar graviton vanishes at leading order the model is stable against linear perturbations Finally, we conclude with a discussion in Section VI. Lastly, for completeness we present the detailed analysis of the Hamiltonian constraint in the Appendix.

II. F (R) NONLINEAR MASSIVE GRAVITY
Following [18], the action of the F (R) nonlinear massive gravity is composed by the UV (F (R) sector) and IR ("massive gravity" sector) modifications, and it writes as where M p = 1/ √ 8πG is the reduced Planck mass, m g the graviton mass parameter, and g µν is the physical metric of the background spacetime. As in the dRGT model, the dimensionless graviton potential U MG is constructed by a group of polynomials through the anti-symmetrization in 4D spacetime as: with In the above polynomials we have introduced the block matrix K ≡ I − X, where X ≡ g −1 f involves a nondynamical (fiducial) metric f . Additionally, [K] denotes the trace of the matrix K. Finally,the coefficients α 3 and α 4 appearing in the graviton potential are the two model parameters.
Similarly to usual F (R) gravity models, it is convenient to perform a conformal transformation on the spacetime coordinate, from the original frame (dubbed as Jordan frame) to the frame where gravity behaves like Einstein's GR (called the Einstein frame). In particular, we impose the conformal transformation as with Ω 2 = F ,R = exp 2 3 ϕ Mp . Correspondingly, the "F (R)" sector of the gravitational Lagrangian can be reformulated as where a potential for the scalar field is effectively introduced: with the comma-subscript denoting differentiation with respect to the following variable.
From the above formulation one could at first sight see a close relation between the F (R) massive gravity and the quasi-dilaton massive gravity [19,20]. However, in the case of quasi-dilaton massive gravity [19] the coefficient in front of the kinetic term of the scalar field is a free parameter, while in our model it is set to unity (this feature has a crucial effect on the perturbational analysis as we will see). Moreover, the quasi-dilaton massive gravity is constructed upon a shift symmetry along the scalar, while this symmetry is in general broken in the present model due to the appearance of the effective potential U (ϕ). Thus, despite the similarity to the quasi-dilaton massive gravity, the two models behave radically differently.
Additionally, the above conformal transformation acts on the gravitational potential sector. Rewriting the graviton potential (2) more conveniently as with and we deduce that in order to perform the conformal transformation (4) we need to incorporate the transformation of the matrix and finally we obtain the transformation as From the expression (12) one can easily see that the gravitational potential in the Einstein frame acquires a scalar-field dependence. This feature is similar to the model of varying-mass massive gravity [21,22]. However, there exist two differences. Firstly, in varying-mass massive gravity the graviton mass is an arbitrary function of the scalar field, while in the present model the scalar-field dependence is fixed by the conformal factor Ω. Secondly and more importantly, in the present model the power index of the scalar-field-dependent (conformal) factor Ω is determined by the order of the gravitational polynomial, while in varying-mass massive gravity it is a common overall factor [21]. In other words, in the present model the separate gravitational terms acquire a different scalar-field dependence, which make the scenario radically different than that of varying-mass massive gravity.
Eventually, to assemble the effects of both the "F (R)" and the "massive gravity" sectors, we can rewrite the Lagrangian of the F (R) nonlinear massive gravity as: The scenario at hand is free of BD ghosts, as can be shown through a detailed Hamiltonian constraint analysis, which for convenience is presented in the Appendix.

III. COSMOLOGY OF F (R) NONLINEAR MASSIVE GRAVITY
In this section we study the cosmological implications of the model under investigation. We consider the standard homogeneous and isotropic Friedmann-Robertson-Walker (FRW) ansatz for the physical metric, while for the fiducial one we start from a Minkowski form f ρσ = η ρσ . Furthermore, we incorporate the usual matter content of the universe, that is we consider an action S m corresponding to an ideal fluid characterized by energy density ρ m and pressure p m respectively, which couples minimally to the gravitational sector. Finally, we fix a particular gauge for the Stückelberg fields.

A. Flat universe
We consider the physical metric in the Jordan frame to be flat FRW: where we have omitted the subscript "g" to distinguish the physical metric, since the fiducial metric is just the Minkowski one. Additionally, we use part of the gauge freedom to impose the following form on the Stückelberg fields: with a c a constant. Inserting the above into the total action (1) we obtain with where the prime denotes the derivatives with respect to t.
For convenience, and in order to present the equations in a compact way, we introduce the following notation: where H is the usual Hubble parameter describing the expanding rate of the universe. Varying the action (16) with respect to all field variables, we obtain the following equations of motion: where we have set N = 1 at the end. In the above expressions we have defined the effective energy density and pressure of the "massive gravity" sector as as well as the effective energy density and pressure contributed by the "F (R)" sector as Here we would like to clarify that, if the model under consideration is used to describe the cosmological dynamics at late times, the effective dark energy component is attributed to the above two contributions, namely Solving the first equation of (19) leads to the following possibilities: H = 0, that is a static universe with a(t) = const; or X 1 = 0, that is A = const, which is a static universe too (a(t) = const). As a result, we find that the F (R) nonlinear massive gravity shares the disadvantage of many other nonlinear massive gravity scenarios, namely that it does not accept flat homogeneous and isotropic solutions (one can easily achieve the same conclusion in the Einstein frame). Consequently, we have to extend the cosmological investigation of this model into the non-flat geometry.

B. Open Universe
Since there exists no acceptable flat solution, in order to study cosmological scenarios we consider an open universe. For completeness, we analyze the cosmological dynamics in both Jordan and Einstein frames.

Analysis in the Jordan frame
We start with the open FRW metric of the form and a k = |K| is associated with the spatial curvature. The Stückelberg scalars take the form of the Milne coordinates: We insert the above formulae into the total action (1) and we obtain where primes denote derivatives with respect to t as introduced in the previous subsection. Similar to the analysis in the flat universe we introduce a series of notation for convenience, namely Then varying the action (29) with respect to the field variables yields the following equations of motion: In the above expressions the effective energy density and pressure of the "F (R)" sector remain the same as in the flat geometry, which are given by (22) and (23), while the effective energy density and pressure of the "massive gravity" sector now become If this model is used to describe late-time cosmology, the effective dark energy component is attributed to the above two contributions through ρ DE ≡ ρ MG + ρ FR and Let us now examine the cosmological equations (32)- (34). Similar to all massive gravity scenarios, Eq. (32) constrains the dynamics significantly. As we observe it leads to two possible solutions. The first one is trivial: a = a k , which leads to the dynamics a(t) = a k t + const. The second one requires Y 1 = 0, and is of physical interest. As pointed out in the self-accelerating backgrounds of the dRGT construction [16], the nontrivial constraint Y 1 = 0 yields This relation can be always fulfilled by choosing b(t) ∝ a(t). In this case expressions (35), (36) imply that ρ MG = −p MG = const, as it is expected similarly to the original nonlinear massive gravity model. This result incorporates the interesting cosmological implication that the graviton mass can induce an effective cosmological constant given by (35), of which the energy density takes the form However, the crucial issue is that in the present model the remaining "F (R)" sector can be chosen at will in the Friedmann equations (33)- (34), leading to a huge class of cosmologies.

Analysis in the Einstein frame
For completeness, we present the investigation of the open geometry in the Einstein frame too. We mention that this will be helpful in the examination of the perturbations of the scenario, since the analysis of cosmological perturbations takes a familiar form in the Einstein frame. We follow the procedure of section II. Specifically, we perform the conformal transformation g µν →g µν = Ω 2 g µν of the open FRW metric (39), with Ω 2 = exp 2 3 ϕ Mp . Therefore, the metric in the Einstein frame is given by and the resulting action in the Einstein frame takes the form withL m = Ω −4 L m describing any additional matter and where primes denote derivatives with respect tot. Note that, as we mentioned in section II, the conformal transformation acts only on the physical metric and not on the fiducial one, and therefore the Stückelberg fields remain unaffected.
We introduce a series of notation for simplicity, namelỹ The equations of motion in this case are obtained by where we have setÑ = 1 in the final expressions, and in the last equation used the relation [3] ∂L As usual, in order to express the results back to the initial metric we use In the above expressions, the effective energy density and pressure of the regular matter component are given bỹ and those of the "massive gravity" sector are expressed asρ Additionally, there exists one more contribution to the total energy density and pressure, namely that of the scalar field, which take the usual forms Finally, in the analysis of the Einstein frame of F (R) nonlinear massive gravity, the corresponding dark energy component is attributed to the "massive gravity" and the scalar field sectors, namelỹ Now let us examine the cosmological equations (43)-(46). Similarly to the analysis of the Jordan frame, Eq. (43) leads to two possible solutions. The first is the trivial solutionH = a k a +Ω Ω , which using (48) is found to be the Einstein-frame equivalent of the trivial solutionȧ = a k found in the Jordan frame. The second isỸ 1 = 0, which further yields This relation can be always fulfilled by choosing b(t) ∝ a(t). Then expressions (50), (51) imply thatρ MG = −p MG , as we found in the Jordan frame too. Note that in the Einstein frameρ MG andp MG depend ont (through the conformal factor Ω), but if we re-express them back into the Jordan frame, in terms of t they become constant. Hence, one can easily deduce that this interesting time varying cosmological constant in the Einstein frame, is due to the conformal factor brought by the coordinate transformation.

IV. PHENOMENOLOGICAL IMPLICATIONS
Having presented the basic background equations of motion of the scenario of F (R) nonlinear massive gravity, we are now able to investigate its phenomenological implications. We mention that, the scenario at hand exhibits a wide class of cosmological behaviors due to the combination of the "F (R)" and the "massive gravity" sectors 1 .
As we discussed in the previous section, the cosmological equations in an open FRW universe governed by F (R) nonlinear massive gravity are given by (32)- (34). The first of these equations constrains the dynamics significantly, leading to a constant contribution of the graviton mass sector, that is ρ MG± = const, given by (38). However, the crucial issue is that in the model at hand, and contrary to usual massive gravity, the remaining "F (R)" sector can be chosen at will in the rest Friedmann equations (33)- (34), leading to very rich cosmological dynamics.

A. The Starobinsky-ΛCDM-like cosmology
One of the most important and interesting sub-classes is when the "F (R)" sector is important at early times and thus responsible for inflation, while the massive graviton is dominant at late times and can drive the universe acceleration as observed today. To give a representative example, we particularly consider the well-known Starobinsky's R 2 -inflation scenario [5] with ξ being the high-order coupling coefficient. This model proves to be the best-fitted scenario with the recently released Planck data [6]. We perform a numerical elaboration of such a cosmological system and in Figure  1 we present the early-time inflationary solutions for four parameter choices. Additionally, in Figure 2 we numerically depict four late-time self-accelerating solutions. We would like to mention that in order to handle the late-time evolutions, we need to define the dark energy equation of state w DE and its density parameter Ω DE . In the F (R)-literature there are two ways to do it. The first is to use the Friedmann equations in the form of (33), (34), in which we introduce an effective Planck mass square M 2 p F ,R , and define Ω FR ≡ ρ FR /(3M 2 p F ,R H 2 ), with w FR ≡ p FR /ρ FR [2]. The second method is to rewrite the Friedmann equations in order to have the standard Planck mass square M 2 p in the left hand side, to absorb the remaining terms in a modifiedρ FR andp FR , and to define  5 (green-doted). All parameters are in units where the present Hubble parameter is H0 = 1, and we have imposed Ωm0 ≈ 0.31, ΩDE0 ≈ 0.69, Ω k0 ≈ 0.01 at the present scale factor a0 = 1 [1]. [3]. Despite the fact that in usual cases the difference is small at the background level (and the specific choice has to be made in accordance with the specific measurement method of M p , Ω DE and w DE ), we stress that the usual conservation equationρ + 3H(ρ + p) = 0 holds only for ρ FR and not for ρ FR [3]. Although, this is not crucial (since we are dealing with modified gravity the dark energy sector is effective, arising from the extra gravitational features, and there is no fundamental reason that it should be conserved), in the present work we follow the second way in order for our model to present energy conservation of the "F (R)" sector, and subsequently for the whole dark-energy one. Thus, the matter, curvature and massive-gravity density parameters are defined using 3M 2 p H 2 too, and they are also conserved. As we observe from these figures, in this particular choice of the Starobinsky ansatz, the Ricci scalar is relative large at early times, and thus the R 2 correction term drives a successful inflation. On the other hand, R becomes very small at late times and thus the "F (R)" contribution is dramatically suppressed, and therefore the acceleration is driven solely by the effective cosmological constant induced by the graviton mass. Therefore, this model provides a unified description of both inflation and late-time acceleration of the universe. In particular, the corresponding cosmological evolution at the background level provides a specific realization of the Starobinsky-ΛCDM cosmology. Note however that at the perturbation level the present model will behave differently comparing to Starobinsky-ΛCDM cosmology, since the graviton sector contributes as well, contrary to the simple cosmological constant which does not. Therefore, strictly speaking, it is a Starobinsky-ΛCDM-like scenario.

B. Late-time Cosmology
As we mentioned above, in the case where the universe evolves at late times, both the "F (R)" and "massive gravity" sectors may contribute effectively to the dark energy component. In particular, the "F (R)" contribution brings dynamical features into dark energy physics, and thus it can be of great phenomenological interest. Thus, in this subsection, we focus on the details of late-time evolutions, going beyond the simple effective cosmological constant of the previous subsection, and using viable F (R) forms, which quantitatively describe the observational data.
We would like to examine the two well-known viable F (R) models that are in best agreement with observations. First we consider the F (R) gravity of the exponential ansatz [24][25][26] with R S and β the model parameters. This ansatz is viable for β > 1 and R S > 0 [24], and it is able to fit observations with only one more parameter comparing to the ΛCDM cosmology. We numerically solve the cosmological system governed by the F (R) massive gravity with this exponential form and we present two different results in Figures 3 and 4, under two different groups of parameter choices, respectively. Our results focus on the dynamics of the equation-of-state parameter of dark energy w DE and the corresponding density parameter Ω DE . and Ω k0 ≈ 0.01 at present [1]. Lower graph: The corresponding evolution of the dark energy density parameter ΩDE (back-solid), as well as its two constituents ΩF R (red-dashed) and ΩMG (green-dotted), of the F (R) and massive-gravity sectors respectively.
For clarity, we also depict the evolutions of both dark energy constituents, namely the separate F (R) density parameter Ω FR and the massive gravity one Ω MG , where Ω DE = Ω FR +Ω MG . For convenience, we use the redshift z = a0 a − 1 as the variable of the horizontal axis, with the present scale factor a 0 set to 1.
In the case of Figure 3 the dark energy is mainly induced by the "F (R)" sector, since the "massive gravity" contribution is quite small. On the other hand, in the case of Figure 4, both the "F (R)" and "massive gravity" sectors have a significant contribution to Ω DE . One notices that in these specific examples, the equation of state parameter w DE evolves from w > −1 into the w < −1 regime and hence provides a realization of the quintom scenario. This is known as a capability and advantage of the F (R) gravity models [27].
Concerning the phantom-divide crossing realization (see [28] for a comprehensive review), there is a wellknown "No-Go" theorem which states that any dynamical dark energy with only one single degree of freedom (DoF) in the frame of General Relativity cannot exhibit it, due to severe gradient instabilities in the dispersion relation at the classical level [29]. In general, in order to circumvent this problems, the quintom models usually involve two scalar DoFs [30], or they are formulated in alternative frameworks such is the non-perturbative stringinspired construction [31], or using a spinor field [32]. Interestingly, in the present case of F (R) nonlinear massive gravity, although at linear level we have only one scalar DoF as we will show in section V, its interaction with the graviton potential modifies the dispersion relation which becomes well-behaved when the phantom-divide crossing occurs. Thus, the aforementioned "No-Go" theorem is avoided (this mechanism is similar to the "No-Go" theorem avoidance in the Horndeski models [33,34]).
The second viable and well-studied F (R) ansatz is of the power-law-like Starobinsky form [35,36]: with λ, n, and R c the (positive) model parameters (n > 0.9 according to Solar-System constraints [36]). We numerically solve the cosmological equations under the above F (R) form and we provide two representative results in Figures 5 and 6, respectively. In the case shown in Figure 5, the dynamics of the dark energy is dominated by the contribution of the "F (R)" sector, since the "massive gravity" contribution is very small. From the evolution of w DE in this figure, one can see that the quintom scenario occurs twice, thus the universe remains in the phantom regime for a short time. This interesting phenomenological property can be applied in the study of very early universe and a nonsingular bouncing solution can be obtained, similarly to [37]. On the other hand, in the case shown in Figure 6, the "massive gravity" sector has the largest contribution at late times, comparing to the F (R) part. Note that the domination of the "massive gravity" sector in Figure 6 makes w DE to approach −1 but does not allow for a crossing. This is expected since, as we mentioned after (37), the pure "massive gravity" sector behaves like a cosmological and Ω k0 ≈ 0.01 at present [1]. Lower graph: The corresponding evolution of the dark energy density parameter ΩDE (backsolid), as well as its two constituents ΩF R (red-dashed) and ΩMG (green-dotted), of the F (R) and massive gravity sectors, respectively.
constant (in the case of asymptotic domination of Ω MG on Ω FR , the w DE tends asymptotically to −1).

V. COSMOLOGICAL PERTURBATIONS
The scenario of F (R) nonlinear massive gravity is free of the BD ghost, as it is shown in detail in the Appendix. However, it is necessary to examine whether this model remain free of instabilities at the perturbative level. The stability issue of various nonlinear massive gravity models has recently been extensively studied in the literature [17,[38][39][40], and the parameter spaces of these models are tightly constrained.
In this section, we proceed to a detailed study of the cosmological perturbations of the scenario at hand, focusing our interest at the level of linear perturbations, and this analysis serves as a completeness of our accompanied work [18]. We work in the unitary gauge which leaves the Stückelberg fields unperturbed. Additionally, as usual we perform the calculations in the Einstein frame, but we drop the tildes on the metric for convenience (thus in this section the dots denote derivatives with respect to the time-variable of the Einstein-frame metric).
The metric variables found in the action, and introduced in Eq. (39), are perturbed according to the following convention The scalar field is also perturbed as ϕ → ϕ + δϕ while the Stückelberg fields are not as we work in the unitary gauge. Therefore, note that in this section the symbol ϕ is used for the scalar field in the Einstein frame, while φ is used for the scalar metric perturbation (thus no confusion is possible since the Stückelberg fields do not appear). Then we can further decompose them in a traceless and transverse tensor (γ ij ), transverse vectors (B i , F i ) and scalars (φ, B, ψ, E) parts through the following relations, where we have introduced the spatial covariant derivative D i compatible with the background 3-metric γ K ij . One can derive the equations of motion of cosmological perturbation by perturbing the extended Einstein's equation of nonlinear massive gravity. We do not include the contribution from extra matter as it does not change the conclusion. Einstein's equation can be expressed as where G µν is the regular Einstein tensor and the extra term X µν comes from the "massive gravity" part of the action. For convenience we introduce which will be used in the following calculation. The process of deriving the perturbed Einstein equations is straightforward. Here we list the non-vanishing equations of motion as follows, which correspond to the (00), (0i), (ii), (ij) components, respectively. Additionally, we have the perturbed Klein-Gordon equation for the field fluctuation, which takes the form of where the coefficientB is given by (42).
Since neither φ norψ appear as coefficient of the graviton mass term, the kinetic structure is very similar to that of GR. This is seen in many other massive gravity constructions, when one considers the self-accelerating cosmological solution as the background. With this in mind, we are allowed to use the Bardeen potentials, to rewrite the equations and reduce the apparent propagating DoFs in the present model. Making use of the above varialbe, Eq. (67) yields Similar to usual perturbation analysis, it is convenient to define the well-known Mukhanov-Sasaki variable: We first combine the (00) (68). In terms of the variable Q, we obtain the perturbation equation as follows: Observing Eq. (72), one deduces that in the scenario at hand the field variable Q is one propagating DoF characterizing the dynamics of cosmological perturbations. The crucial issue is to examine whether ψ B is a dynamical DoF too, or whether it can be integrated out through some constraint equation.
We first re-express equations (64) and (65) in terms of the Bardeen potentials (69) and Q. Therefore, we can derive the following two useful formulae: with the operator C defined as Equations (73) and (74) determine the dynamics of ψ B . However, they are not independent since d dt ψ B =ψ B . The consistency of their combination, together with the main perturbation equation (72), leads to a strong constraint, namely E = 0. Hence, we can clearly see that there exists only a single scalar propagating DoF, that is Q, at the linear order in the perturbation theory around the cosmological background.
Let us now examine the stability of the single DoF Q. Inserting (73) into (72), we observe that the sign in front ofQ remains equal to one, and thus no apparent ghost instability appears. The existence of the term containing ψ B would modify the dispersion relation of the gradient terms. But in a realistic cosmological background the spatial curvature term is observationally suppressed and thus the corresponding effect on modifying the dispersion relation would be secondary. The last interesting property shown in Equation (72) is the term proportional to m 2 g . Obviously, this is a new term brought by the gravitational potential and thus it is interpreted as the effective mass of the scalar graviton. Despite the corrections from other terms, one deduces that ifỸ Q < 0 then the model is free of tachyonic instabilities. Note that this requirement still leaves a large part of the parameter space, and thus a huge class of interesting cosmological behavior.
We close this section by making some comments on the comparison with other extended nonlinear massive gravity models, such are the quasi-dilaton massive gravity [19] and the mass-varying massive gravity [21], following the discussion of section II. From the point of view of perturbation analysis, both the quasi-dilaton and massvarying massive gravity models involve at least two scalar DoF, one introduced by the additional scalar field and the other being the longitudinal mode of the graviton. This can be verified by counting the number of nonzero eigenvalues of the matrix for the kinetic part of the perturbation action [39]. Applying the method of [39] in the present scenario of F (R) nonlinear massive gravity (this can be performed immediately by taking the results of [39] and set the coefficient in front of the scalar-field kinetic term to unity), we find that there exists only one nonzero eigenvalue in our model, and this implies only a single DoF. Thus, the present model is radically different from the aforementioned two.

VI. CONCLUSIONS
Recently, a nonlinear massive gravity theory was constructed in [11], in which the BD ghost can be removed in the decoupling limit to all orders in perturbation theory, through a systematic construction of a covariant nonlinear action. The theoretical advantages of such a construction led to a wide investigation of its cosmological implications [15,16], as well as of the black hole solutions [41], while the connections to bi-metric and multi-metric gravity were also revealed [42,43]. However, the cosmological perturbations around these solutions are found to exhibit in general severe instabilities [17], and thus extensions of the basic theory are necessary.
In this work, we constructed the F (R) extension of nonlinear massive gravity, as an accompanied paper to our previous work [18], focusing on a detailed cosmological investigation, both at the background and perturbation levels. In particular, we investigated three representative cosmological evolutions. In the first one we considered the F (R) sector to be of the Starobinsky's R 2 -form and we explicitly showed that this model can describe both inflation and late time acceleration in a unified picture, with the F (R) sector driving inflation and the graviton mass responsible for dark energy. Furthermore, we studied the case in which F (R) takes the usual viable exponential form and we showed that this class of models can easily quantitatively describe the late-time universe behavior, including the possibility of a phantom-divide crossing. In the last example we considered the powerlaw-like Starobinsky F (R) ansantz, which is another viable form, with similar quantitative behavior. In these two classes of evolutions, the effective dark energy component includes contributions from both sectors, namely the F (R) and the massive graviton one.
Apart from the interesting cosmological implications at the background level, the advantage of the present scenario is the behavior of cosmological perturbations. We developed the perturbations by expanding the Einstein equations to linear order, around a cosmological solution. Although the dispersion relation of the scalar mode is different from the case of GR, the effect is negligible in a realistic cosmological background due to the suppression of the curvature term. Our analysis showed that there exists only one propagating scalar mode in this model, which can additionally be free of ghost instabilities in a large part of the parameter space. This issue is the crucial difference of F (R) nonlinear massive gravity comparing to other extensions, such are the quasi-dilaton and the varying-mass massive gravity. The fact that both F (R) and massive graviton sectors are needed for this behavior, may imply that the UV and IR modifications of gravity may not be independent.
At the end of the present work, we would like to highlight one interesting property of the specific model of R 2massive gravity. One notices that, the effective potential for the scalar field after conformal transformation, approaches a constant in the UV regime. This extremely flat potential, which can be applied to drive sufficiently long inflationary process at early universe as analyzed in Sec. IV A, indicates an approximately shift symmetry along the scalar field. This feature interestingly coincides with the model of quasi-dilaton massive gravity with a specifically chosen dilatonic parameter. However, when the scalar field evolves into the IR regime, it is stabilized at the vacuum and hence this shift symmetry can be spontaneously broken. This profound property is interesting to the process of symmetry breaking in nonlinear massive gravity models, and deserves further investigations.
Amongst all the coefficients in the metrics, the lapse N 0 g and the shift N g (the three N i g 's expressed as vector) of the physical metric, as well as the corresponding ones for the fiducial metric, N 0 f and N f respectively, are not dynamical DoFs.
In general massive gravity γ ij allows for at most six propagating modes, one of them being the origin of the BD ghost. Thus, a potentially healthy theory must maintain a single constraint onγ (from now on a bar denotes the matrix form) and the corresponding conjugate momenta, along with an associated secondary constraint, which will lead to the elimination of the ghost DoF. In the following we show the existence of these constraints in F (R) massive gravity.
In general, the propagating DoFs are the {γ ij , ϕ}, such that their phase space spans 14 dimensions after including their conjugate momentas {π ij , π}. In standard GR, four ghostly propagating DoF are eliminated by the constraints from the equations of motion of the lapse N and shifts N i . The inclusion of a generic mass term changes these equations and they do not directly give rise to any constraints on γ ij . However, the dRGT mass term is such that the four equations of motion for N and N i give rise to an equation independent of N and N i , which provides a constraint that eliminates one potentially dangerous DoF. Consistency of the first constraint to be valid at later times introduces a secondary constraint. These are necessary since we expect massive gravity to have only five propagating DoF carried by γ ij . The extra 6-th mode carried by the helicity 0 mode is at the origin of the BD ghost.
We can make this constraint explicit by introducing a Lagrange multiplier in the action. This is done by using a new shift function n i defined through the following relation: where D i j is determined by The explicit form of D is not needed for the analysis, but without loss of generality we apply the expression that satisfies the identity ξD = Dξ. Defining λ = (1 − n i ξ ij n j ), the Hamiltonian can be written as with A =β 1 λ 1/2 +β 2 λD i i + n i ξ ij D j k n k +β 4 ξ γ One can observe that the Hamiltonian (A3) is not linear in n i and hence the equations of motion δL δn i = C i = 0 do not lead to any constraints on the propagating DoFs, but they determine n i (γ ij , π ij , ϕ, π). However, (A3) is now linear in N , and thus δL N = C = 0 provides our first constraint on the system.
Consistency requires that C = 0 is valid at all times. If this is trivially satisfied, then there would be no addi-tional constraints. Fortunately, we find that the derivative does not vanish trivially, and this indicates that it does imposes the secondary constraint needed to eliminate the ghost. In the above we have introduced the Poisson brackets First of all, we need to prove that {C(x), C(y)} = 0. Otherwise, the appearance of the lapse N in (A5) would not turn this consistency requirement into a constraint on the propagating modes as: C (2 ∂γmn (x). The involved poisson brackets read In order to deal with the delta functions, the smoothing functions f (x) and g(x) are introduced as follows,