Viable cosmology in bimetric theory

We study cosmological perturbations in bimetric theory with two fluids each of which is coupled to one of the two metrics. Focusing on a healthy branch of background solutions, we clarify the stability of the cosmological perturbations. For this purpose, we extend the condition for the absence of the so-called Higuchi ghost, and show that the condition is guaranteed to be satisfied on the healthy branch. We also calculate the squared propagation speeds of perturbations and derive the conditions for the absence of the gradient instability. To avoid the gradient instability, we find that the model parameters are weakly constrained.


I. INTRODUCTION
Our universe consists of various particles with different spins and masses. The graviton, the spin-2 particle mediating the gravitational force, is of special interest since gravity is the least understood among fundamental forces in nature. Assuming Lorentz invariance, Weinberg's theorem in 1964 [1] and its extensions [2,3] exclude more than one interacting massless gravitons in four-dimensional Minkowski spacetime. Those theorems, however, do not exclude massive graviton(s) interacting with a massless graviton.
In the present paper we consider a bimetric theory of gravity, i.e. a physical setup involving two dynamical metrics interacting with each other. In this setup, after diagonalizing the mass matrix for metric perturbations around Minkowski background, we end up with a massive graviton and a massless graviton, in accord with the above mentioned general theorems. Bimetric theory thus propagates seven physical degrees of freedom in Minkowski background: five from the massive graviton and two from the massless graviton. Until recently, however, it was thought that nonlinear extension of massive gravity inevitably would have involved a sixth degree of freedom (eighth degree of freedom in nonlinear bimetric theory, e.g. [4]), which would have been a ghost [5]. This ghost degree of freedom, called Boulware-Deser (BD) ghost, was recently excised in the construction of a massive gravity theory by de Rham, Gabadadze and Tolley (dRGT) [6,7] at fully nonlinear level [8,9]. A simple extension of the dRGT massive gravity allows the construction of a fully nonlinear bimetric theory of gravity without the (would-be) BD ghost [10]. It is thus this formulation that the studies of bimetric theory in the present paper are based on.
Having a promising candidate for theoretically consistent bimetric theory, it is important to investigate whether it can accommodate viable cosmology. Before starting the study of cosmology in bimetric theory, however, let us briefly review the current status of cosmology in dRGT massive gravity and its extensions.
In the covariant formulation of dRGT massive gravity, the basic quantities in the gravity sector are a metric field and four scalar fields called Stückelberg fields. The original dRGT theory respects the Poincaré symmetry in the space of Stückelberg fields so that the Stückelberg fields enter the action only through the so-called fiducial metric, which is the pullback of Minkowski metric in the field space to the spacetime. It turned out that the dRGT theory with the internal Poincaré symmetry does not allow for any non-trivial flat Friedmann-Lemaître-Robertson-Walker (FLRW) solutions [11]. The same no-go result holds for closed FLRW solutions as well. One can nonetheless find non-trivial and self-accelerating open FLRW solutions in this theory [12]. Slightly extending the theory by replacing the Minkowski fiducial metric with de Sitter or FLRW one, it also becomes possible to find not only open but also flat and closed FLRW solutions with [13] or without [14][15][16] self-acceleration. Unfortunately, all those FLRW solutions in dRGT theory turned out to be unstable due to either linear instability called Higuchi ghost [17] or a recently found nonlinear instability [18,19].
The origin of the new nonlinear instability found in [18] is the fact that kinetic terms of three among five degrees of freedom in FLRW backgrounds are exactly proportional to the equation of motion for the temporal Stückelberg field and thus vanish on shell [13]. Those kinetic terms vanish at the quadratic level in the action, but do appear at the nonlinear level and can become either positive or negative, depending on the nature of perturbations.
One can in principle evade the instability of cosmological solutions in dRGT massive gravity by relaxing the FLRW symmetry, i.e. either homogeneity [11] (see also [20][21][22][23][24][25] for related solutions) or isotropy [19,26]. Another possibility is to extend the theory by introducing an extra scalar field in the gravity sector [27][28][29]. In both cases, the above mentioned exact proportionality between the kinetic terms and the background equation of motion is detuned and thus the nonlinear instability can in principle be avoided. The nonlinear bimetric theory of gravity can be considered as yet another extension of dRGT massive gravity, due to an extra spin-2 field. Hence, the above mentioned no-go result for stable FLRW cosmological solution does not directly apply to the nonlinear bimetric theory. Unfortunately, in nonlinear bimetric theory the analogue of the self-accelerating branch still suffers from nonlinear instability. We thus study the stability of FLRW solutions in the other branch in the bimetric theory originally proposed in the context of the graviton oscillations and gravitational wave observations [30], which we call the healthy branch. 1 In this study, we take into account two matter fields each of which couples to either the first or second metric.
The rest of the present paper is organized as follows. In Sec. II we describe the model of our interest. We briefly review the background cosmology assuming the spatial homogeneity and isotropy in Sec. III, where we also introduce the notion of two different branches. Then, we discuss linear perturbations around this cosmological background in Sec. IV. Starting with the pure gravity case, we discuss the tensor, vector and scalar-type perturbations one by one. We identify the conditions for the absence of ghost and gradient instability. Section V is devoted to summary and discussion.

II. MODEL
The covariant action for the gravity sector is constructed out of two four-dimensional metrics g µν and f µν . It is the sum of the two Einstein-Hilbert actions I EH,g and I EH,f for the metrics g µν and f µν , respectively, and the non-derivative mixing term I mix which is built by requiring that the Boulware-Deser ghost is absent at all orders [8]. Including the matter, the total action we consider is where Each term in I mix is constructed as where the square brackets denote trace operation and The square root in this expression represents the matrix that satisfies and whose eigenvalues are all positive. Notice that the cosmological constant term for f µν can be expressed as a linear combination of L i (i = 0, . . . , 4), while that for g µν is L 0 . Besides the three mass scales M g , M f and Max i |α i |m, we thus have four dimensionless parameters (four among α i (i = 0, . . . , 4)). This means that, for each choice of the set of mass scales M g , M f and Max i |α i |m, we have freedom to tune four additional quantities, e.g. ξ c , J(ξ c ), J ′ (ξ c ) and Λ/m 2 defined in the next section. As for the matter sector, we consider scalar fields φ g and φ f minimally coupled to g µν and f µν , respectively. Because of their equivalence to irrotational barotropic perfect fluids (see e.g. [34]), we restrict our consideration to k-essence fields 2 , with the action with When we simply write P a , with a = g, f , it denotes the background value of P a (X a ) and corresponds to the background pressure in the perfect fluid picture. We also introduce the energy density and sound speed of the effective fluid as where prime denotes derivative with respect to the argument.

III. BACKGROUND EQUATIONS AND SOLUTION BRANCHES
In this section, we derive the equations of motion for the background FLRW universe, and also discuss the branches of the background solutions. We denote the two metrics g µν and f µν as with where N = N (t) and n = n(t) are the background lapse functions, and a = a(t) and α = α(t) are background scale factors. Notice that the two background metrics are both homogeneous and isotropic, having common isometries. The background equations are given by 3 1 where we have defined H ≡ȧ/(N a), H f ≡α/(nα), κ ≡ M 2 f /M 2 g and an overdot represents derivative with respect to t. Here, we have also introducedρ where ξ andc are, respectively, the ratios of the background scale factors and lapse functions defined as and For later convenience, we define In the next section we shall see that m eff is the effective graviton mass. 3 Combining Eqs. (13), (15) and (17), or equivalently, Eqs. (14), (16) and (18), we obtain a constraint, This constraint indicates that there are two branches of solutions: one is specified by the condition J = 0 and the other by H = ξH f . Before starting the analysis of the healthy branch, we briefly mention the other branch with J = 0. Since the condition J = 0 implies that ξ = α/a is constant, the quantitiesρ m,g andρ m,f are also constant. In this case the quadratic order of perturbation of the mixing term of the action, I mix , is equivalent to the cosmological constant terms on both metrics. (See Eq. (45) below with J being set to 0.) Thus, the Hamiltonian structure for linear perturbations is the same as that in two copies of general relativity. As a result, instead of the 7 degrees of freedom we expect from a healthy bimetric theory, we end up with 4 degrees of freedom that are dynamical at linear order. [36]. This branch is an analogue of the self-accelerating branch in massive gravity [12], which is known to possess incurable problems, such as a less-than-expected number of the degrees of freedom at quadratic level [13] and appearance of a non-perturbative ghost [18]. Thus, we consider only the healthy branch specified by H = ξH f in the rest of this paper.
First, we derive two algebraic relations, i.e. constraints that hold in the H = ξH f branch. Combining Eqs. (13) and (14) under the assumption H = ξH f , we have where we have introduced a function of ξ defined bŷ Equation (25) should be interpreted as the equation that determines the value of ξ. The other constraint is derived from d(H − ξH f )/dt = 0. Using Eqs. (15) and (16), this constraint can be rewritten as where we have defined This constraint is to be interpreted as the equation that determines the difference between the light cones of two metrics,c − 1. Hence, barring special tuning of model parameters, vanishing of (ρ g + P g ) − ξ 2 (ρ f + P f )/κ means the crossing ofc = 1 and W keeps a definite sign along the trajectory of the time evolution. The regimec > 1 is preferred as a viable model when the ordinary matter fields are coupled to the g metric. This is because, wheñ c < 1, the electro-magnetic wave travels faster than the propagation speed of the f -metric perturbations and hence a UHECR traveling with a speed very close the speed of light may emit the f -gravi-Cherenkov radiation, which is severely constrained by observations [37,38]. In Ref. [30] a healthy background cosmology was proposed. We generalize the analysis in [30] to the case with two matter fields, removing also the restriction to the low energy regime. In this solution at low energies, ξ converges to a constant ξ c (assumed to be ξ c = O(1)) that solveŝ The mass scale of the coupling term m is assumed to be parametrically large, but the effective graviton mass m eff (ξ c ) is kept small, by tuning the model parameters.
The other constraint (27) impliesc → 1 at low energies. Namely, the difference of the light cone between two metrics vanishes in the low energy limit. We expandρ m (ξ) around ξ = ξ c aŝ where we defined Hence, Eq. (25) implies that Using this, we find that the Friedmann equation (13) can be approximated as where we definedM Eq. (34) can be interpreted as the Friedmann equation for two matter fields ρ g andκ −1 ρ f with the effective gravitational constantM 2 g and the effective cosmological constant Λ. We did not assume smallness of |Λ/m 2 | so far. The pure gravity case can be easily obtained by taking the limit of ρ a = 0 and P a = 0. In this case both metrics are de Sitter and we have ξ = ξ c andc = 1.
If we tune the effective cosmological constant Λ so that Hereafter, we assume Eq. (36) as well as Eq. (29) to hold when we take the low energy limit. We remark that we do not intend to solve the cosmological constant problem in the present paper and that the condition (36) can be realized by simply tuning the α 0 parameter. In the low energy limit or in the pure gravity case, W in Eq. (27) with K = 0 is reduced to m 2 eff /2 − H 2 , where m eff is the effective graviton mass defined in Eq. (23). This quantity must be positive for the absence of Higuchi ghost [17]. As mentioned above, the sign of W does not change in general. Therefore, we choose the branch in which is satisfied. In this case, the condition to avoid the Cherenkov radiation,c > 1, is reduced to The positivity of W also indicates that J = 0 is not realized except in the limit where H 2 +K/a 2 and W simultaneously vanish, provided that |K|/a 2 < H 2 in accord with observation. However, in this low energy limit, the value of ξ converges to ξ c , which is different from the zeros of J(ξ) in general. In other words, J(ξ c ) = 0 unless fine-tuned. Hence, we find that the healthy branch does not cross the J = 0 branch. The positivity of W may break down only when the sequence of solution disappears as we increase the energy scale. From the constraint (27), we find thatc − 1 diverges when W crosses 0. Notice that we can also expressc − 1 as Hence, whenc − 1 diverges,ξ also diverges. This indicates that the system would exit the regime of validity of the effective field theory there. By contrast, when the right hand side of Eq. (27) vanishes, it just meansc − 1 crosses zero in general. Therefore, this does not mean the flip of the sign of W . Hence, we conclude that, in the regime of validity of the effective field theory, the condition W > 0 is maintained. The above phenomena can be understood more intuitively when there is no matter field coupled to the f -metric. Notice that W is related to dρ m /d ln ξ as Hence, when there is no matter coupled to the f -metric, W > 0 is identical to dρ m /d ln ξ > 0, which is the condition for the absence of Higuchi ghost discussed in Ref. [39]. The meaning of this condition is clear: as we increase the matter energy density ρ g , the constraint (25) with ρ f = 0 implies thatρ m should decrease. When the minimum of the functionρ m (ξ) is reached, we cannot extend the background solution beyond that critical energy density of ρ g . Conversely, as long as the healthy branch solution continues to exist, the condition dρ m /d ln ξ > 0 is maintained.

IV. PERTURBATIONS AROUND FLRW BACKGROUNDS
A. Pure gravity case We expand the two metrics perturbed around FLRW backgrounds as where (Φ, V i , H ij ) and (ϕ, v i , h ij ) represent perturbations. For the present discussion, it is convenient to decompose H ij and h ij further into their trace parts and traceless parts as where TrH ≡ γ ij H ij , H T ij ≡ H ij − 1 3 γ ij TrH, Trh ≡ γ ij h ij , and h T ij ≡ h ij − 1 3 γ ij Trh. The Einstein-Hilbert action for g µν is expanded as and the quadratic part is given by where A similar expansion applies to the Einstein-Hilbert action for the f -metric. Up to second order, the mixing term (4) is expanded as In the pure gravity limit, the quadratic action can then be diagonalized by introducing new perturbation variables Using the relations for the background, we can rewrite the quadratic action for perturbation as with where The effective gravitational coupling for the "+" fields, M 2 + , the effective cosmological constant, Λ, and the effective graviton mass for the "−" fields, m eff , have been already defined in Eqs. (37), (32) 4 and (23), while the effective gravitational coupling for the "−" fields is given by In the absence of matter, we find that the linear combination of metric perturbations described by "+" fields is just the linearized general relativity on a de Sitter background, while the other linear combination described by "−" fields forms a decoupled massive spin-2 field with mass m eff around the same background.

B. Inclusion of matter
Having finished the analysis on the pure gravity case, we study the quadratic action taking the matter sector into account. Then, in the subsequent subsections we study the tensor, vector and scalar sectors in turn, and argue the stability conditions for each of them.
We introduce the perturbation of the matter fields as We expand the action for the matter fields (9) up to second order as and the second order term is given by The action for the φ f field can be obtained similarly.

C. Tensor sector
We start with our analysis on the tensor sector. We restrict the perturbation of the three metrics to the transversetraceless perturbations as To keep the notation simple, we suppress the superscript T T below. Combining all terms (two Einstein Hilbert terms, matter terms and the interaction terms), the action quadratic in tensor modes reduces to We stress that up to now, we did not assume any branch and only used the background equations (15) and (16). In this most general setup, H ij and h ij fields have a generically time-dependent coupling. Taking the low energy limit where ξ ≃ ξ c andc ≃ 1, the tensor action can be put into the diagonal form where H ± ij are defined as in Eqs. (47). This action is essentially the same as the tensor part of Eq. (49) for the pure gravity case. H + ij is the massless graviton mode, and H − ij is the massive graviton mode with mass m eff given by Eq. (23).

D. Vector sector
We introduce vector perturbations to the metric through where B i , b i , E i and S i are transverse with respect to the covariant differentiation associated with γ ij metric. Using the background equations (13)-(16), the total quadratic action of the vector perturbations becomes In the above form, the action is manifestly gauge invariant, since Next, we vary the action with respect to the non-dynamical degrees B i and b i , to obtain At this stage, it is convenient to expand the modes in vector harmonics. The solutions to the above equations then read where we defined and k 2 is the eigenvalue of −∆. Substituting the solution for the auxiliary modes (61) into the action (59), we obtain where we defined the gauge invariant combination and the squared propagation speed of the vector mode One can immediately verify that A = 0 on the J = 0 branch and thus the vector modes are non-dynamical at linear order. For the healthy branch, A should be positive at k → ∞ to avoid the ghost instability, and this condition is satisfied if J > 0. This latter condition is automatically satisfied since we have chosen the branch that satisfies (38).
If we take the low energy limit, the squared propagation speed of perturbation, c 2 V , is reduced to As discussed in Ref. [30], d ln J/d ln ξ is necessarily large for the Vainshtein mechanism to efficiently work in this model. However,c − 1 is suppressed at low energies. In fact, when Γ is dominated by the first term of the right hand side of (22), i.e., when Γ ≈ ξJ, the absolute value of the right hand side of the above equation is always less than unity and c 2 V is thus guaranteed to be positive. On the other hand, if we consider the case with |d ln J/d ln ξ| > 1, i.e., if the right hand side of (22) is dominated by the second term, then the right hand side of the above expression is O(m 2 eff /m 2 ξJ), which can be larger than O(1). If this is the case, the model parameters are constrained to satisfy d ln J/d ln ξ > 0 to avoid the gradient instability.
In the main text, we did not take into account the rotational modes of fluids. In Appendix, the case with the matter coupled to the g-metric only is discussed taking into account the rotational modes.

Reduction of the quadratic action
The scalar perturbations are introduced through in the metrics, and in the matter sector. We integrate out the non-dynamical modes B, b, Φ and ϕ, ending up with an action depending on six variables, ψ, Σ, E, S, δφ g and δφ f . As we have not fixed the gauge yet, there are two pure gauge degrees in this action. Furthermore, we also expect that the would-be Boulware-Deser mode should be non-dynamical. Hence three of the six variables are to be eliminated. Under the coordinate transformation where ξ µ g = (ξ 0 g , ∇ i ξ g ), the six variables transform as where The three physical degrees of freedom can be made manifest by choosing the following gauge invariant variables: If we adopt the so-called flat gauge conditions E = ψ = 0, which fix the gauge completely as seen from Eq. (70), Y 1 and Y 3 coincide with the modes δφ g and S, respectively. In this gauge, after expressing δφ f in terms of Y 2 , the Boulware-Deser mode is manifestly non-dynamical in the action and can be integrated out. The resulting action after expanding the fields into harmonics is where Y = (Y 1 , Y 2 , Y 3 ) T is the field array, while K T = K, N T = −N and M T = M are 3 × 3 real matrices 6 .

No-ghost conditions
We first discuss the conditions for avoiding ghost. The kinetic matrix can be diagonalized, by applying the rotation 6 In this paper we do not discuss the effective Newton potential, partly because linear analysis is not sufficient to study the metric perturbation within the Vainshtein radius. Here we quote the result for the Newton potential within the quasi-static analysis in the linear perturbation, when the matter is coupled only to the g-metric. The resulting Newton potential δΦ is given by as Then, we find that the ghost is absent if the following inequalities are satisfied whose explicit expressions are where To avoid the catastrophic ghost instability, these conditions must be satisfied, at least, in the k → ∞ limit. In this limit conditions (77) reduce to These conditions are satisfied if The first two conditions are just the null energy conditions for the matter fields, which are requested for the sound waves of fluids to be stable. Let us now focus on the last condition, B > 0. We can easily verify that B can be rewritten as Therefore, as long as the branch with W > 0 is concerned, the positivity of B is guaranteed. In the low energy limit or in the pure gravity case, we find where m eff is the effective graviton mass defined in Eq. (23). Thus, one can see that the condition B > 0 is the extension of the Higuchi bound for the absence of ghost modes.

Sound speeds
We now study the speeds of propagation of the scalar modes, in the high frequency limit. The equation of motion is obtained from the action (73) as Assuming monochromatic time-dependence of perturbation Y ∝ e −i ωN dt and the adiabaticity of the backgrounḋ ω/N ≪ ω 2 , we obtain the dispersion relation The eigenfrequencies can be found by solving this equation. Expanding the eigenfrequencies around k → ∞ as ω 2 ∼ c 2 s k 2 /a 2 , we can read the coefficients c 2 s , the squared propagation speeds of perturbation in the high frequency limit, as These must be positive to avoid a fatal instability. The positivity of the first one is not obvious. It would be instructive to rewrite it as This expression proves that c 2 s = 1 in the pure gravity case. In the low energy limit, the leading order terms come from the first line on the right hand side of the above equation and they reduce to Then, the positivity of c 2 s gives a constraint on the model parameters in a similar sense to c 2 V . Before ending this section, we compare our result with the solutions discussed in [36]. In this case,ρ m is dominated by the term proportional to 1/ξ, so here (and only here) we assume ξ ≪ 1. As a result, we have J ∼ −κξρ m = constant, hence d ln J/d ln ξ ≈ 0. If we assume that there is no matter field coupled to the f -metric, and hence the gradient instability is inevitable. This instability occurs in the regime of validity of the low energy effective field theory and thus is physical for m ≪ H ≪ Λ 9 (≪ Λ 3 ), where Λ n = (M g m n−1 ) 1/n . 7 If a matter sector which couples to the f -metric is present and is the dominant contribution of the right hand side in Eq. (25), a similar calculation yields Again, for both non-relativistic and radiation fluids, the scalar sector suffers from gradient instability.

V. SUMMARY AND DISCUSSION
We have presented a linear analysis of cosmological perturbations in bimetric theory, in which two metrics are coupled through non-derivative coupling so that it does not yield Boulware-Deser ghost as prescribed in [10]. We consider perturbations around the background of two dynamical FLRW metrics sharing spatial isometries, each of which is minimally coupled to a different k-essence field. The contracted Bianchi identity of either metric yields a constraint which defines two branches. The J = 0 branch, in which the ratio of the scale factors of two metrics is fixed to a constant value determined by the condition J = 0, contains only four dynamical degrees of freedom at linear order [36], while at non-linear order, is known to suffer from instabilities [18]. In this work, we focused on the other, at least, seemingly healthy branch, and studied the stability of the cosmological background against perturbations.
In the absence of matter fields, the linearized action is found to be a combination of two decoupled spin-2 fields: one of which is linearized GR, while the other is a linearized spin-2 theory with a mass term. In the setup where two matter fields are minimally coupled to respective metrics, we considered a specific scenario where the mass scale of the interaction term between two metrics is large. When we choose this mass scale to be small of the order of the expansion rate H, Refs. [32,36] concluded that either crossing a singularity or encountering an instability is inevitable. We adopted the basic idea proposed in [30], where the singularity and the instability are expected to be outside the reach of the low energy effective theory since the mass parameter is chosen to be large enough. Even in such a setup, the extra forces can still be screened under certain conditions, owing to the Vainshtein mechanism. In this paper we showed that we can actually avoid the singularity and the instability, at least, at the level of linear perturbation around FLRW background.
For the general setup, we identified the necessary conditions for avoiding ghost instabilities. We also derived the expression for the propagation speed of each perturbation mode and the condition for the avoidance of gradient instabilities. We found that the matter perturbations have positive kinetic energies as long as the background matter fields do not violate the null energy condition, and that they propagate at the ordinary sound speeds of the corresponding fluids. From the positivity of the kinetic term of the gravitational degree of freedom in the scalar sector, we obtained the extension of the Higuchi bound [17] (see also [40]), which turns out to be automatically satisfied as far as the healthy branch is concerned. We also found that the squared propagation speeds of gravitational degrees of freedom in the scalar and vector sectors can be negative, which leads to the so-called gradient instability. In order to avoid this instability, the model parameters are weakly constrained.
The way how matter couples to gravity is always an issue in alternative gravity theories. In the present paper we have considered two matter fields and supposed that each of them minimally couples to one of the two metrics. This seems a particularly safe choice: in the massless limit m → 0, the system is decomposed into totally decoupled two subsystems, each of which consists of a massless graviton interacting with a matter field. This is trivially consistent with Weinberg's theorem in 1964 [1] and its extensions [2,3] that exclude more than one interacting massless gravitons, since the two subsystems are completely decoupled from each other in the massless limit. While other possibilities for matter coupling remain unexplored, any theoretically consistent schemes of matter coupling should be in accord with the above mentioned general theorems. For example, Ref. [41] proposes that a matter field can couple to both metrics. In the massless limit of this prescription it appears that the two massless gravitons can interact with each other through the matter field. It is interesting to see whether and how this can be reconciled with the general theorems.
Finally, the present study should help generalizing the predictions of the bimetric theory on future experiments (see e.g. [43]).
Note added: Shortly after this work was completed, Ref. [44] appeared on the arXiv, where solutions with effective mass of the order of the present-day Hubble rate are considered (see also [45][46][47][48]). We remark that our analysis does not exclude healthy backgrounds other than the ones considered in the present work. a 10 m 2 JξM 2 g (ρ + p) 16N 4 (c + 1) > 0 .
On the chosen background these conditions imply J > 0. This condition is automatically satisfied by the healthy branch background solution. We find two different speeds of propagation. Namely, one is c 2 V in Eq. (65) and the other is The second mode corresponds to the degree of freedom in the matter sector.