On the viability of bigravity cosmology

We revisit the question of viability of bigravity cosmology as a candidate for dark energy. In the context of the low energy limit model, where matter couples to a single metric, we study linear perturbations around homogeneous and isotropic backgrounds to derive the Poisson's equation for the Newtonian potential. Extending to second order perturbations, we identify the Vainshtein radius below which non-linear scalar self interactions conspire to reproduce GR on local scales. We combine all of these results to determine the parameter space that allows a late time de-Sitter attractor compatible with observations and a successful Vainsthein mechanism. We find that the requirement on having a successful Vainsthein mechanism is not compatible with the existence of cosmological solutions at early times.


I. INTRODUCTION
Einstein's theory of general relativity (GR) [1] has been the widely accepted theory of gravity with the impeccable ability to match observations for over a century [2]. However, the discovery of the accelerated expansion of the Universe [3] has lead to the construction of many modified theories of gravity which attempt to account for this observation in a more natural way than the addition of a cosmological constant to the Einstein field equations. See Refs. [4][5][6] for reviews.
Dropping the notion of a massless spin-2 graviton is arguably the most natural extension of GR. The effect of endowing the graviton with a nonzero mass was first considered in 1939 by Fierz and Pauli in a linear construction [7]. In this theory, the mass term is built by requiring the absence of negative energy states (ghosts) and breaks the linearized diffeomorphism invariance, thus resulting in a massive spin-2 field theory propagating 5 degrees of freedom (d.o.f.). Following this initial work, van Dam and Veltman [8] and, independently, Zakharov [9] showed that in the massless limit of the linear theory, GR is not recovered. The cause of this discrepancy is that the helicity-0 component of the graviton does not decouple from the trace of the stress tensor of the matter source in this limit. The resolution to this problem, offered by Vainshtein, is to extend the theory to include nonlinear self-interactions to allow a smooth GR limit [10].
However, nonlinear theories generically suffer from the Boulware-Deser ghost, the unwanted 6th d.o.f. which breaks the linear tuning of Fierz and Pauli [11]. The ghost-free potential was constructed after almost 40 years, by de Rham, Gabadadze and Tolley (dRGT) in a decoupling limit [12] and was subsequently shown to be ghost free to all orders [13]. The theory is built out of a tensor ffiffiffiffiffiffiffiffiffi ffi g −1 f p constructed from a physical metric g μν and a flat fiducial metric f μν . In the context of cosmology, massive gravity with dRGT potentials does not allow exact cosmological solutions without generating pathologies. The self-accelerating branch suffers from nonlinear ghost instability [14], while the normal branch does not allow expansion [15] in the case of a flat fiducial metric. Extending the flat fiducial metric to maximally symmetric space-times [16], the normal branch can support cosmology. However, in the case of de Sitter fiducial metric [17], these either suffer from a Higuchi ghost [18] or do not have a successful Vainshtein mechanism [19], while for anti-de Sitter, the cosmology is protected against acceleration [20].
The situation is less severe for the bigravity theory, where the f μν metric is promoted to be dynamical [21]. In the model in which matter is coupled only to the physical metric g μν , although the self-accelerating branch cannot evade the conclusions of Refs. [14], the normal branch can sustain cosmology. For a late-time acceleration that is sourced by a massive graviton, the mass needs to be generically of the order of Hubble rate today. For this scenario, Ref. [22] showed that the scalar perturbations in the radiation dominated era suffer from a gradient instability, ruling out a viable cosmology. 1 There have been various studies on the stability of this model [25][26][27][28][29], and two ways to circumvent this conclusion have been proposed. The first is to impose a hierarchy between scales, effectively decoupling the massive graviton from the matter sector, thus making the model indistinguishable from GR [30]. The second way is the so-called "low energy limit" [31], where the bare mass parameter is allowed to be large m ≫ H 0 , while the late-time accelerated expansion can be achieved via a fine-tuning of coupling constants. This tuning introduces a hierarchy between the interaction parameter m and the effective mass of the dynamical graviton modes. The stability of this model was shown for a large portion of the parameter space in Ref. [32], while some implications for primordial gravitational waves were studied in Ref. [33].
In this work, we focus on the low energy limit model in which matter is coupled to a single metric. Although the model has interesting gravitational-wave phenomenology [31], its implications for cosmology have not been explored in detail. The goal of the present work is to fill this gap and determine the late-time implications of this model. The paper is organised as follows. In Sec. II, we review the bigravity theory and give the equations of motion. In Sec. III, we discuss the background evolution, focusing on the late-time de Sitter attractor. In Sec. IV, we study both linear and second order scalar perturbations, to obtain the analogue Poisson's equation and to identify the Vainshtein radius, respectively. In Sec. V, we summarize all the conditions obtained to fix three model parameters and discuss potential observable signatures in this reduced parameter space. We conclude with a discussion in Sec. VI.

II. BIGRAVITY THEORY WITH DRGT INTERACTIONS
In this section, we give a brief review of bigravity. The action for the theory is given by [21] S ¼ M 2 where α n are dimensionless parameters and M g and M f are the corresponding Planck scales for the two metrics g and f. In the above, L n are the dRGT interaction terms given by L 0 ðKÞ ¼ 1; where square brackets denote trace operation and we defined the building block tensor as We note that the square root above is a tensor operation defined by the relation In this setup, g μν corresponds to the physical metric, i.e., the metric that matter sector L matter couples to, while f μν is a dynamical background field. We now vary the action (1) with respect to g μν and f μν , to yield the equations of motion for the g and f metrics, respectively, where G μν and G μν are the Einstein tensors built out of the g and f metrics, respectively, and is the energy-momentum tensor for the matter sector. We show the explicit result for the variation of L n with respect to g μν and f μν in Appendix A. Using (A3), one can verify that (5) matches the expressions in Ref. [34].

III. BACKGROUND COSMOLOGY
In this section, we study the background cosmology under the Ansatz that both metrics take Friedmann-Lemaître-Robertson-Walker form in the same coordinate system, where nðtÞ and NðtÞ denote the lapse functions, while aðtÞ and αðtÞ represent the scale factors for the g and f metrics, respectively. For the matter content, we consider a perfect fluid described by the energy-momentum tensor, where u μ is the 4-velocity of the fluid and satisfies the condition u μ u μ ¼ −1, P is the pressure, and ρ is the energy density. In accordance with the homogeneous and isotropic metric Ansätze, the background values for the pressure and energy density are functions of time only, while In what follows, we will restrict our discussion to a matter sector consisting only of a pressureless nonrelativistic fluid with P ¼ 0.
For the background metrics (7) and a perfect fluid (8), the equations of motion (5) reduce to four independent equations [32], where a dot represents a time derivative and we defined the following functions in accordance with the notation of Ref. [32], with 4 , and we also have From here on, we set α 2 ¼ 1 without loss of generality, effectively absorbing this parameter into the definition of m. The contracted Bianchi identity for individual metrics yields an effective constraint. For instance, differentiating (9) then combining it with (11) and (13), we obtain This equation branches out into two solutions, where the self-accelerating branch J ¼ 0 is known to lead to nonlinear ghost instabilities [14]. Instead, we choose the normal branch with H ¼ ξH f . This solution links the evolution of the f-metric to the g-metric one, and the consistency of the two Friedmann equations gives an algebraic relation between ξ and the matter density ρ, where we defined the combination In order to avoid the early time gradient instability in this branch [29], we will adopt the low energy limit defined by which allows us to push the instability beyond the reach of the effective field theory [31,32]. As time evolves, the density ρ redshifts as a −3 , and the solution for ξ converges to a constant value ξ c defined bŷ To describe the evolution close to this late-time attractor, we linearize Eq. (17) around ξ ¼ ξ c to relate the departure from this point to the matter density, where Λ ≡ m 2 ρ m;g ðξ c Þ and the subscript c corresponds to the values of functions evaluated at the de Sitter attractor. Following Ref. [32], we now assume jΛ=ðm 2 J c Þj ≪ 1.
Using these results, we find that the Friedmann equation (9) can be approximated as VIABILITY OF BIGRAVITY COSMOLOGY … PHYS. REV. D 99, 104032 (2019) where the effective Planck scale isM 2 g ≡ ð1 þ κξ 2 c ÞM 2 g , and we now identify Λ as the effective cosmological constant. 2

IV. COSMOLOGICAL PERTURBATIONS
In this section, we consider perturbations around the low energy background model and determine the effect of the two-metric interaction on the linear growth of structure. We outline the method taken to study the perturbations in theory, the process to isolate the scalar mode in the Poisson's equation, and the extension to nonlinear order. We only consider scalar perturbations in this work as they are the only relevant ones for the large scale structure.

A. Linear perturbations
To begin, we perturb both metrics around the backgrounds (7) in the Newtonian gauge for the g metric, where ðϕ; ψ; b; S; ϕ f ; ψ f Þ are the perturbation variables and we fixed the time coordinate such that N ¼ 1. The perturbations in the matter sector are introduced via ρ ¼ ρðtÞ þ δρ and u μ ¼ ð1 − ϕ; ∂ i vÞ, giving With these decompositions and using the quasistatic approximation [35], we can derive an analogue of Poisson's equation for the potential ϕ [32], where k corresponds the momentum of the mode in the plane-wave expansion. The derivation is summarized in Appendix B. The contribution from the two-metric interaction is encoded in the function W defined by This quantity also plays a major role in the perturbative stability conditions, with W > 0 corresponding to the bigravity generalization of the Higuchi bound [32]. The form of Poisson's equation is similar to the one in the presence of a scalar field source. The traceless part of the g equations of motion that it is given by (B3), reveals that the perturbation S acts as a source for anisotropic stress.

B. Vainshtein radius
We now move on to the study of second order perturbations and an identify the scale at which the perturbative expansion breaks down. This will allow us to determine the Vainshtein radius where the scalar graviton decouples from the matter sector and the evolution closely follows GR.
In order to do this, a few approximations are in order. In addition to restricting the study to the de Sitter attractor, we consider scales where the expansion can be neglected. We also focus on small scales and assume ∇ 2 ≫ m 2 . Keeping metric perturbations up to second order in the equations of motion, we find that, unlike the g-metric equations, the nonlinear f-metric equations do not exhibit the enhancement ∇ 2 =m 2 with respect to the linear part. This allows us to solve the linear f-metric equations and substitute the solutions for ðψ f ; ϕ f ; ϕÞ into the nonlinear components of the g-metric equations. The solutions are Upon substitution of Eqs. (28) into the g-metric equations, we obtain where we kept only the second order terms that are enhanced in the limit ∇ 2 =J c m 2 ≫ 1. Using Eq. (29) to replace ∇ 2 ψ, Eq. (29) reduces to where the nonlinear term has the expected Galileon-like structure. The scale at which the nonlinear terms become important depends on the normalization of the field. We defineS such that the linear traceless equation of motion (27) reduces toS ¼ ψ þ ϕ. With this normalization, the nonlinear equation (31) becomes Thus, the nonlinear term dominates for C∇ 2S ∼ Oð1Þ, revealing the order of the Vainshtein radius as where r g corresponds to the Schwarzschild radius of a spherical body. This result is consistent with the similar calculation in Ref. [31].

V. FIXING MODEL PARAMETERS
We are now ready to fix the model parameters without introducing the bare cosmological constant term, which is equivalent to setting α 0 ¼ 0. On the other hand, we keep α 1 nonzero. This means that the theory does not admit the Minkowski solution, which is not a problem as we are interested in cosmological solutions.
We start by trading α 3 for ξ c using the definition of the de Sitter fixed point (20) and obtain α 3 ¼ α 3 ðξ c ; α 1 ; α 4 ; κÞ as We then fix α 4 matching the effective cosmological constant Λ in the approximate Friedmann equation (22) to the observed value, which is equivalent to solving The solution for α 4 ¼ α 4 ðξ c ; α 1 ; κÞ is Finally, the last parameter is fixed by requiring a sensible Vainshtein radius which ensures that the effects of modified gravity are hidden below a certain distance scale to recover GR on solar system and galactic scales. The workings of the Vainshtein mechanism are that derivative selfinteractions of the scalar are enhanced around a matter source such as a star, so that the effect of the fifth force is screened below the Vainshtein radius. The relation we will consider is We introduce the parametrization where b ∼ Oð0.1 − 1Þ [36,37] and allows us to tune the size of the Vainshtein radius. The equation which relates the model parameters to the Vainshtein radius is (35) From this relation, we then fix α 1 , where we made use of the fact that H 0 ≪ m. Noting that the solution with the þ sign leads to a negative W [defined in Eq. (26)], we will choose the − sign solution in the following. Moreover, the solution only exists for ξ c > 1.
After this procedure, we reduce the number of free parameters down to two: ξ c and κ. We now check whether there are any inconsistencies in the background equations of motion. Expanding the left-hand side of Eq. (17) around the attractor, we have In the limit H 0 ≪ m, the coefficients of the linear and quadratic terms are In Sec. III, we used the linear term to obtain the approximate Friedmann equation (22). However, we see that the first derivative is suppressed by H 0 =m, while the second derivative term is manifestly positive. As a result, when the quadratic term dominates, there is no real solution to this equation. This observation allows us to determine the parameter range, which grants a physical evolution. The linear term is dominant if in which case the solution to (43) behaves as Using the condition (45), the above relation yields an upper bound on b; Since the matter density today is of order of Therefore, the solution exists for Turning this relation around, given a parameter b, the cosmological description can go as far back as a in , before which no physical evolution exists. Although we set α 0 ¼ 0 in order not to introduce a bare cosmological constant, we can check that this conclusion holds even if α 0 ≠ 0.
Suppose we wish to describe the evolution of the scale factor up from the last scattering surface onward. Therefore, we set a in ¼ a CMB ¼ 10 −3 . In order to have the low energy limit (19) be valid at the time of the cosmic microwave background (CMB), the minimum mass parameter allowed is m ¼ H CMB , where H CMB is the Hubble parameter at the last scattering surface. From Eq. (49), this initial value of the scale factor corresponds to a value of b ¼ 10 −9=2 . In order to check this estimate, we compared the exact numerical solution of Eq. (17) to the linear approximation. The comparison is summarized in Fig. 1. The exact solution only appears around a ∼ 10 −3 , after which the value of ξ becomes closer to the de Sitter attractor value ξ c .
We close this section with a discussion of the effect on the large scale structure. From Eq. (25). we can determine the consequences of the several tunings; in order to have an observable effect, the quantity W has to be comparable to the k 2 contribution. The function W can be interpreted as the effective mass of the gravity perturbations and behaves as W ∼ m 2 J. It encodes the information about the scale at which the modifications to gravity appear. Using the approximate expression, we can estimate its value, using m ∼ H CMB ∼ 10 9=2 H 0 , which implies that the effect of the two-metric interaction appears only at scales smaller than approximately 0.1 Mpc where linear perturbation theory is no longer applicable.

VI. CONCLUSIONS
We have presented an analysis of the linear and nonlinear perturbations in bigravity where nonderivative two-metric coupling is introduced as in Ref. [21] so as not to generate the Boulware-Deser ghost. We considered a perfect fluid with equation of state P ¼ 0 coupled to the g metric and studied metric perturbations around Friedmann-Lemaître-Robertson-Walker, while adopting the healthy branch of solution with H ¼ ξH f . Poisson's equation was derived at the linear level in perturbations, and we identified the modification to Poisson's equation due to the extra d.o.f. present from the massive graviton. Furthermore, we studied perturbations going beyond linear order and identified the Vainshtein radius, below which the derivative self-interactions of the scalar screen the effect of the fifth force and conspire to reproduce GR on local scales.
We then looked at the effect that fixing three of the model parameters has on the background cosmology of the theory using the following requirements: ensuring the existence of the late-time de Sitter attractor, matching the effective cosmological constant in the Friedmann equation to the value we observe, and proposing we have a Vainshtein radius. Bigravity in the low energy limit can admit a sensible cosmological solution, but this comes with the cost of lowering the Vainsthein radius. With a value of b ¼ 10 −9=2 , for which the Vainshtein radius is given by R 3 v ¼ 10 −9 H −2 0 r g , we are able to describe the evolution of the scale factor up until the last scattering surface at a in ¼ 10 −3 . However, to satisfy the observational constraints, we need to impose b ¼ Oð0.1-1Þ, which results in a very short window of the viable cosmological evolution.
The main conclusion of this paper is that the stable bigravity model that is distinguishable from GR does not provide a reasonable description for the late-time acceleration of the Universe. With this result, we have established that none of the exact cosmological solutions to dRGT massive gravity/bigravity theory, where matter couples to a single metric, admits a viable and testable dark energy model. There are several potential directions from here. The most conservative option is to explore extensions which preserve the local Lorentz invariance of dRGT, while avoiding new dynamical d.o.f. One such extension is to allow matter to couple to a composite metric, which avoids the generation of Boulware-Deser ghost within the range of the effective field theory [38]. In the normal branch of cosmological evolution [39], the vector modes suffer from a gradient instability in the radiation era [40], while the scalar mode becomes a ghost in the late Universe [41]. The selfaccelerating branch, which would be problematic in the noncomposite theory, becomes detuned by the effect of matter and allows a stable evolution that undergoes a bounce [41,42]. For the composite coupled theories, if the double coupling extends to the Standard Model sector, this implies that the light cone corresponding to the observed gravitational waves from binary neutron star merger [43] is different than the one for photons [44].
Another extension that persists is the generalized massive gravity, where the translation symmetry in the Stückelberg scalar field space is broken, while the Lorentz invariance remains intact [45]. In this construction, the number of d.o.f. is the same as in dRGT, although the parameters are promoted to functions of the scalar fields. This theory admits self-accelerating open universe solutions, and their stability was shown in the decoupling limit. However, its phenomenology remains largely unexplored.  [46] appeared and claimed that our conclusions on the viability of bigravity cosmology cannot hold in general. As we clearly stated in the abstract and the rest of the paper, we are considering the low energy limit of the theory, and our aim was never to find a general conclusion; nor did we make this claim. Both of the counterexamples that were highlighted in Ref. [46] (i.e., hierarchy between the Planck scales and screening of gradient instabilities) are already mentioned in the Introduction.

APPENDIX A: DERIVATION OF EQUATIONS OF MOTION
In this Appendix, we compute the variation of the interaction term with respect to g μν and f μν . We define Using this definition, we can vary the trace of various powers of this tensor, δ½X n ¼ n 2 ðX n Þ α μ ðg αν δg μν − f αν δf μν Þ; which is valid for any power n ≥ 1. In the above, we made use of δf μν ¼ −f μα f νβ δf αβ . The variation of the interaction terms can be written in the following form: X α μ ðg αν δg μν − f αν δf μν Þ; ðg αν δg μν − f αν δf μν Þ; ðg αν δg μν − f αν δf μν Þ;