Bigravity and Lorentz-violating Massive Gravity

Bigravity is a natural arena where a non-linear theory of massive gravity can be formulated. If the interaction between the metrics $f$ and $g$ is non-derivative, spherically symmetric exact solutions can be found. At large distances from the origin, these are generically Lorentz-breaking bi-flat solutions (provided that the corresponding vacuum energies are adjusted appropriately). The spectrum of linearized perturbations around such backgrounds contains a massless as well as a massive graviton, with {\em two} physical polarizations each. There are no propagating vectors or scalars, and the theory is ghost free (as happens with certain massive gravities with explicit breaking of Lorentz invariance). At the linearized level, corrections to GR are proportional to the square of the graviton mass, and so there is no vDVZ discontinuity. Surprisingly, the solution of linear theory for a static spherically symmetric source does {\em not} agree with the linearization of any of the known exact solutions. The latter coincide with the standard Schwarzschild-(A)dS solutions of General Relativity, with no corrections at all. Another interesting class of solutions is obtained where $f$ and $g$ are proportional to each other. The case of bi-de Sitter solutions is analyzed in some detail.


Introduction
In recent years, many proposals to modify gravity at very large, i.e. cosmological, distances have been put forward with the hope to find an alternative explanation to the observed acceleration of the Universe. Indeed, one can imagine that this acceleration is not due to some, yet unknown and dark, component of the Universe, but to some failure of General Relativity at large distances (see e.g. [1]).
To obtain such a modification of gravity, one of the simplest ideas -in fact already present in the original motivation of Einstein to introduce a cosmological constant -is obviously to give some sort of mass to the graviton. However, it has long been known that theories of massive gravity are rather problematic. In the linearized regime, a Lorentz invariant mass term for the graviton must have the Fierz-Pauli (FP) form [2], or else ghosts appear in the spectrum. With the FP mass term, light deflection caused by static sources differs from that of General Relativity (GR) even in the limit of massless gravitons, a pathology known as the van Dam-Veltam-Zakharov (vDVZ) discontinuity [3].
Among the recently studied proposals to modify gravity at large distances, Dvali-Gabadadze-Porrati gravity ( [4] DGP gravity) has been the subject of numerous works in particular because of its ability to produce a late time acceleration of the expansion of the Universe even with a vanishing cosmological constant [5]. Another theory, sharing with DGP gravity the same tensorial structure of the graviton propagator, and (for that reason) also the vDVZ discontinuity, is non linear massive gravity. The latter theory can be obtained from some bigravity theory, where one of the two metrics is frozen (see [6] and references therein). On the other hand, bigravity theories have the ability to produce an accelerated expansion of the Universe in a non standard way [7], one of the two metrics being then regarded as some new type of dark energy component ( it is interesting to note that there is an integration constant appearing in the effective cosmological constant in these models, which reminds one of unimodular gravity [8]). The vDVZ discontinuity, without a cure, would suffice to rule out any theory in which it appears. It was however proposed that this discontinuity could disappear non perturbatively, i.e. that it would not be seen in exact solutions of the theory, but only shows up as an artifact of the linearization procedure [9]. While this seems problematic in non linear massive gravity [10,11], there is some evidence that it does work in DGP gravity [12,13,14]. It is then of some importance to better understand exact solutions of various types of massive gravities, among which bigravity theories. In this paper we concentrate on spherically solutions of the latter theories. Indeed, in bigravity theories, non-derivative interactions between two different metrics can give a mass to some of the spin-2 polarizations and a large class of non trivial exact solutions are known. This is not the case for more complicated models such as DGP gravity, and thus bigravity theories provides an invaluable framework to investigate various properties of "massive" gravity and related models of large distance modification of gravity.
Interestingly there is another motivation to take a closer look at bigravity theories. Indeed, it has been shown [15,16] that the vDVZ discontinuity can also be avoided when the mass term breaks Lorentz invariance. The different phases of massive gravity with explicit breaking of Lorentz invariance have been further analyzed in [16]. Some examples with spontaneous breaking (due to additional vector field condensates) have also been considered in Ref. [17] (see also [18]). The construction in [17] exhibits ghosts and tachyons at low momenta, but it is nevertheless phenomenologically acceptable for certain choices of parameters. The question of Lorentz violation consistent with cosmological observations has been recently studied in [19]. Such Lorentz violating mass terms also appear in the so-called ghost condensate scenario with non trivial cosmological consequences [20]. It seems worthwhile, in this context, to explore alternative scenarios with spontaneous breaking of Lorentz invariance. A straightforward example is provided again by spherically symmetric solutions of bigravity theories.
In Ref. [21] we studied the global properties of a wide class of bigravity solutions, and here we will develop the theory of linearized perturbations around some of them. A particularly interesting example is the solution where both metrics f and g are flat, but with different values of the "speed of light". This is the simplest case where Lorentz invariance is spontaneously broken. It is also interesting in its own right, since it corresponds to the far field in a wider class of spherically symmetric exact solutions, of the Schwarzschild form.
The paper is organized as follows. In Section 2 we review the basics of bigravity. We show that if the interactions between the two metrics are non-derivative, then there are always a broad class of solutions in the Schwarzschild-(A)dS family, with Lorentz-breaking asymptotics. In Section 3, we study the perturbations to Lorentz breaking bi-flat solutions. In Section 4 we consider perturbations around solutions where the two metrics are proportional to each other, concentrating in the case of bi-deSitter solutions. Section 5 summarizes our conclusions.
While this paper was being prepared, a related work [22] appeared where the study of perturbations to bi-flat Lorentz-breaking backgrounds is outlined. Where we overlap, our findings agree with those of Ref. [22] (see also [23] and [24]).

Exact Solutions
Following [25], we consider the action Here L f and L g denote generic matter Lagrangians coupled to the metrics f and g respectively, and subindices f and g on the Ricci scalar R indicate which metric we use to compute it. For the background solutions, we shall restrict attention to the case where there is only a vacuum energy term in each matter sector L f = −ρ f , L g = −ρ g , where ρ f and ρ g are constant. The kinetic terms are invariant under independent diffeomorphisms of the metrics f and g, but the interaction term is invariant under "diagonal" diffeomorphisms 6 , under which both metrics transform. There is much freedom in the choice of the interaction term in (1). For instance, Ref.
[25] considered a non-linear generalization of the Fierz-Pauli mass term with As we shall see in the following Section, linearized perturbations around asymptotically bi-flat Lorentz-breaking solutions of this particular theory show a singular behaviour. For that reason, it will be instructive to look at generalizations of (2). The most general interaction potential which preserves the "diagonal" diffeomorphism takes the form [6] where τ n = tr[M n ], n : 1, ..., 4 correspond to the traces of the first four powers of M µ ν = f µα g αν , and V is an arbitrary function.
There is also some arbitrariness in the way one introduces matter fields, since one has two different metrics at hand. One can at least couple matter minimally to any one of the two metrics, in which case the other metric can be regarded as some kind of exotic new type of matter. Those two choices correspond to the two matter Lagrangians L g and L f , of action (1), where it is understood that the matter fields entering into L g and L f are different. This opens the possibility to have two types of matter, one which feels the metric g and the other which feels the metric f . In fact one can imagine more complicated situations in which matter fields would be coupled to some composite metric built out of the two metrics f and g. If one wishes to recover the standard equivalence principle, one should obviously ask that standard matter only couples to one metric, and a minimal choice is thus, e.g., that all matter fields appear say in L f (respectively L g ), while L g (respectively L f ), will be simply given by a cosmological constant. With such a choice, matter moves along geodesics of the metric f (respectively g), and, provided the solutions for the metric f are the same as in standard general relativity (which turns out to be possible as will be seen below), there would be no deviations from 4D general relativity seen in matter motion. However, various question arise, should one wish to consider bigravity theories as realistic. Among those, the fact that bigravity theories are likely to suffer from the same type of instability (at the non linear level) discovered by Boulware and Deser in non linear massive gravity [26,27] (see however [6,7]). Our viewpoint is here more to use those theories as a toy to study properties of more complicated models, than pushing the idea that bigravity theories are realistic. Thus we will not further discuss in the following phenomenological consequences of various types of couplings to matter. In the aim of discussing the vDVZ discontinuity, we will only consider the possibility to have two types of matter coupled respectively minimally to metrics f and g, as indicated in the action (1).
Let us then introduce the exact solutions subject of this study. The general static spherically symmetric ansatz for bi-gravity can be written as [28] g µν dx µ dx ν = Jdt 2 − Kdr 2 − r 2 dθ 2 + sin 2 θ dφ 2 (5) where the metric coefficients are functions of r. Note that in general it is not possible to write both metrics in diagonal form in the same coordinate system.
For the potential (2), and for D(r) = 0, it was shown in [28] that the general solution is given by where Here β is an arbitrary constant, γ = 2/3 and the potentials p and q are functions of r. Somewhat surprisingly, these potentials turned out to coincide with those of the standard Schwarzschild-(A)dS family. Solutions of the form (7-8) are called Type I.
Substituting the Type I ansatz into the effective energy momentum tensors which are obtained by varying S int , one readily finds that these take the form of cosmological terms: whereΛ Note that the corresponding cosmological constants for the metrics (7-8) are not determined solely by the vacuum energies ρ f and ρ g . They also contain a contribution from the interaction term in the Lagrangian. This contribution depends not only on the parameters ζ and u (recall that v = 1/2 − u), but also on the arbitrary integration constant β, as seen from (12)(13).
A crucial observation [21] is that even without assuming any specific form for the functions p(r) and q(r), the Type I ansatz leads to energy momentum tensors of the form (11). Here, we would like to clarify the reason for that, and to show that Type I solutions (as well as some generalizations thereof) exist also in the generic case (4). Indeed, for arbitrary metrics f and g where we have introduced the notation where l is the number of derivatives. Moving to the frame where both metrics are diagonal, the matrix M = f −1 · g can be put to the diagonal form with eigenvalues λ i . Two arbitrary metrics g µν and f µν which are solutions of the vacuum Einstein equations will be solutions for bi-gravity if all the τ n are constant and the eigenvalues of the matrix entering (14)(15) are all equal to each other. Note that for a given ansatz, the constancy of the traces (or of the eigenvalues) is a frame independent notion. The equations of motion will be then satisfied provided that where Λ f and Λ g are the cosmological constants.
Remarkably, the non-trivial background (7)(8) has the property that the eigenvalues of M are constant Thus, it is enough to impose In the frame where M is diagonal the previous combination is a constant diagonal matrix with only two different constant eigenvalues Both eigenvalues will coincide when This tells us that for any potential there will exist non-trivial solutions with certain γ and β satisfying (19) (note that the values of V (n) depend also on β and γ). In addition, (17-18) must be satisfied. These are three equations for the parameters Λ f , Λ g , β and γ. Therefore, it is clear that one of the effective cosmological constants can be chosen arbirarily. It has the status of an integration constant.
Another interesting class of solutions is obtained by taking f and g proportional to each other, but otherwise arbitrary In this case, the matrix M is proportional to the identity M µ ν = γ −1 δ µ ν and the energy-momentum tensors (14-15) readΛ Thus, for any matter content this term just adds to the vacuum energy. From Bianchi identitiesΛ f andΛ g must be constant, and f and g must then be solutions of the vacuum Einstein's equations. Generically, the expressions forΛ f,g depend on γ, so that they imply a constant γ. In this case, the parameter γ is determined through Einstein's equations by noting that (20) implies Clearly, this class will include solutions in the Schwarzschild-(A)dS family, although non-spherically symmetric solutions are possible as well. Note also that such solutions can easily be generalized to multigravity theories by deconstructing 5D metrics with a warp factor [29].
Let us now turn to the study of perturbations to some of these background solutions.

Perturbations around Lorentz-breaking biflat metrics
In a theory with two metrics with Einstein-Hilbert kinetic terms and no interaction, there are 4+4 ADM Lagrange multipliers. When we add a non-derivative interaction which preserves diagonal diffeomorphisms, 4 combinations of these may in principle appear non-linearly in the action [6]. For these, their equation of motion relates them to the other variables, but they do not lead to further constraints. Thus, we have a minimum of 4 and a maximum of 8 Lagrange multipliers for 20 metric components. Hence, we generically expect a maximum of (10 − 4) + (10 − 8) = 6 + 2 = 8 degrees of freedom and a minimum of (10 − 8) × 2 = 2 + 2 = 4. In a Lorentz-invariant context, the first possibility corresponds to a massless and a massive graviton, whereas the second would correspond to two massless gravitons. In the Lorentz breaking context, it is possible to have a massive graviton with just two physical polarizations [30,31].
Let us consider a general potential V [{τ n }] as in (4). As we showed in the last section, the vacuum energies ρ f and ρ g can be tuned so that the previous potential has asymptotically bi-flat solutions. At large distances from the origin, these take the form whereη and η µν = diag(1, −1, −1, −1). The parameters γ and β are related by Eq. (19). For β = 1 7 , we cannot simultaneously write both metrics in the canonical form η µν , and Lorentz invariance breaks down to spatial rotations. It will be convenient to introduce the general perturbation in the form whereη µν is the inverse ofη µν . The perturbation to the metric f has been defined with the upper indices, just because this simplifies the manipulations which yield the action quadratic in the perturbations shown below. For the remainder of this section, all space-time indices will be raised and lowered with the canonical Minkowski metric η µν . The interaction Lagrangian quadratic in perturbations then reads where, after imposing (19) For the sake of simplicity, we will restrict to potentials V [{τ n }] for which Eq. (19) is independent 8 of β, and determines γ. From equation (28), this implies n 0 = n 4 . In particular, this class includes the interaction (2), which, as we shall see, leads to a rather pathological behaviour for the perturbations. On the other hand, it is general enough to be representative of generic choices of potentials. In Refs. [15,16] the case of a single graviton with a Lorentz violating mass term has been discussed. For comparison with those references, it will be useful to introduce where c > 0 is an irrelevant constant which has the dimensions of mass squared. Note that the components h g 0i and h 0i f are absent from (27). As noted in [22] the absence of such terms is a consequence of invariance under diagonal diffeomorphisms in this background. In the case of a single graviton (with a Fierz-Pauli kinetic term), the absence of h 0i in the mass term leads to a very interesting behaviour [16,30,32], where the two polarizations of the massless graviton acquire mass, while all the other modes do not propagate. 9 Let us now investigate whether a similar phenomenon occurs in our model. The situation is not directly reducible to that of a single graviton, since the equations of motion are not diagonal. Also, the kinetic term breaks the Lorentz invariance. It is convenient to decompose the perturbations into irreducible representations of the spatial rotations, where t X ii = t X ij,i = V X i,i = F X i,i = 0 for X = f, g, and all space-time indices are raised and lowered with the metric η µν .
To second order in the perturbations, the kinetic terms in (1) can be written in terms of these scalar, vector and tensor variables as: where˜ =η µν ∂ µ ∂ ν ,κ f = γ −1 β 1/2 κ f and dot means a derivative with respect to time. At the linear level, the transformations generated by independent diffeomorphisms δx µ = ξ µ X in each one of the metrics can be expressed as Note that the kinetic term is written in terms of the following quantities: which are invariant under both gauge symmetries. On the other hand, the full action (including the mass terms), is invariant only under the diagonal gauge symmetry No second order scalar combination of h X 0i is invariant under this gauge symmetry, which implies that those terms are always absent (cfr. (27)). We may now analyze the propagating degrees of freedom.

Tensor Modes
The linearized Lagrangian for the tensor and vector modes can be expressed as whereκ f = γ −1 β 1/2 κ f . The corresponding equations of motion in Fourier space read from which we obtain the dispersion relations where κ 0 = n 2 (βκ g +κ f ) and κ 1 = n 2 (κ g +κ f ). At high energies, we have In this limit, each one of the two gravitons propagates in its own metric (with the corresponding "speed of light" 10 ) along null directions k µ = (ω, k) satisfying The low energy expansion of (38) is given by The first dispersion relation corresponds to two massless polarizations which propagate at the "intermediate" speed Note that for β > 1 we have β −1 < c 2 s < 1, while for β < 1 we have 1 < c 2 s < β −1 . The second dispersion relation, Eq. (41), corresponds to two massive polarizations. It is easy to check that the graviton polarizations are stable and tachyon free as long as κ 0 > 0, in the whole range of momenta k. The second dispersion relation (41) corresponds to the massive graviton.

Vector Modes
From the Lagrangian (35), we find that V g i and V f i do not appear in the interaction term. Varying with respect to the vector fields we have, We can always use the diagonal diffeomorphism invariance to work in the gauge where V g i = 0. It then follows from (42) where F i are arbitrary functions of position and f i are arbitrary functions of time. The latter are in fact irrelevant, because F X i enters the metric only through spatial derivatives. Formally, we may describe this as a gauge symmetry , which we can use in order to write, without loss of generality, It then follows from (43) that , where again we eliminate the additive time dependent part. Finally, from (44) we obtain wheref i are new arbitrary functions of time. This is not a desirable situation, since it means that the initial conditions do not determine the future evolution of V f i . Technically, the absence of the fields V g i and V f i in the mass term leads an enhanced gauge symmetry in the linearized Lagrangian. Indeed, we can consider independent gauge transformations for each of the metrics of the form ξ X i = ξ X i (t). As we have discussed, these do not affect the F X i , but can be used to give both of the V X i an arbitrary time dependence.

Scalar Modes
The Lagrangian for the scalar modes can be expressed as Let us first study the non-homogeneous modes. The mass terms do not depend on B g nor on B f , so those fields are Lagrange multipliers, just as in Einstein's gravity. Variation with respect to these fields yields ∆ψ g = ∆ψ f = 0.
The variation with respect to A g and A f yields the constraints Once we substitute the first of these constraints in the Lagrangian, the quadratic term in E h and E l takes the form ( We can now distinguish two different cases, neither of them with propagating scalar degrees of freedom. First, if the coefficient n 2 −n 0 +n 3 does not cancel, the equations of motion for E h and E l result in a new constraint which determines these fields, and upon substitution into the Lagrangian we are left without any scalar degrees of freedom. If the coefficient cancels, as happens for the potential (2), E g and E f are Lagrange multipliers appearing in the gauge invariant combination E h + E l . After using (48), the variation with respect to E h yields ∆ψ g = ∆ψ f = 0.
The Lagrangian cancels after substitution of these constraints, and there are no propagating degrees of freedom. Note that in this last case the combination E h + E l , is not determined by the equations of motion. Again, this is not a desirable feature, since it means that the value of this combination, which is gauge invariant under the diagonal diffeomorphisms, is not predicted by the linear theory. Nevertheless, we expect that higher order terms in the expansion will determine E h + E l , since there is no symmetry in the non-linear Lagrangian under which this quantity can be "gauged" to arbitrary spacetime dependence. Concerning the homogeneous modes, after using the constraints we are left with two modes ψ f and ψ g which have a negative definite kinetic term. Nevertheless, the dispersion relations for the degrees of freedom which diagonalize the equations of motion are ω 2 = 0 and ω 2 = M 4 n 2 (κ f +κ g ) > 0, so there is no classical instability associated to these modes.

Coupling to Matter and vDVZ discontinuity
All known explicit and non-singular exact solutions of bigravity (which we reviewed in Section I) are also solutions of GR. This immediately suggests that the vDVZ discontinuity may be absent altogether in this theory at the non-linear level. Also, from the analysis of perturbations done in the previous Section around the Lorentz breaking background, it is clear that the situation here is very different from that of ordinary massive gravity. The massive spin-2 graviton has only two physical polarizations (as opposed to the five polarizations of the ordinary FP massive graviton), and there are no propagating vectors or scalars.
Let us consider the coupling of the linearized theory to conserved sources. To this end, we introduce the couplings where T µν g and T f µν are conserved, i.e. ∂ µ T µν g = 0 and η ρµη ρα ∂ α T µν f = 0. In terms of the decomposition (30), we have where we have introduced the gauge invariant combinations Inverting the equations of motion for the tensor modes in the presence of the source T ij , we find and an analogous expression for t f ij : In the limit M 4 → 0 this reduces to the standard expression for linearized GR. For the vector modes, the equations of motion read ∆ It follows immediately that ∆(F g i + F f i ) = 0, and therefore the term proportional to M 4 vanishes. This means that there is no difference with the GR results for each one of the metrics.
For the scalar part, we may start with variation with respect to B X i , which yields the constraintṡ and where A + = A g + β −1 A f , ψ + = ψ f + ψ g , and E + = E f + E g . Variation with respect to ∆E X yields, with the help of (57), Substituting into (59), we have and using (57), we have For (n 2 − n 0 + n 3 ) = 0, this determinesĖ + in terms of the sources. The solution will depend on an arbitrary time independent mode E 0 (x). For the singular case (n 2 − n 0 + n 3 ) = 0, Eq. (61) do not determine E + at all. Instead, it imposes some non-trivial equations to be satisfied by the sources, which seem hard to motivate. Thus, coupling to the sources seems rather inconsistent in this case, unless (n 2 + 3n 3 − 3n 0 ) = 0 as well. But this would imply n 2 = 0, in which case the tensor modes are massless.
In the generic case, the solution for the ψ potentials is of the form where C + ( x) is entirely determined by initial conditions. Finally, variation with respect to ψ f and ψ g leads [after use of (63)] to the following equations for the gauge invariant potentials: where (66) In general, the solution depends on an arbitrary "initial" function C + ( x). This corresponds to a mode with dispersion relation ω 2 = 0 in the linear theory. It was argued in [16] that in such cases, from higher order terms the expected dispersion relation will be of the form ω 2 ∼ p 4 , and in this sense C + corresponds to a slowly varying "ghost condensate" [34]. In what follows, we shall take the initial condition C + ( x) = 0.
For n 2 − n 0 + n 3 = 0, the solution is of the form and Hence, there is a well behaved massless limit, with corrections of order M 4 ∆ −2 to the gauge invariant potentials Φ and ψ. This means, in particular, that there is no vDVZ discontinuity. This is quite analogous to the "half massive gravity" model discussed by Gabadadze and Grisa [31] (see also [30]). The additional terms lead to corrections to the Newtonian potential. The sign of this correction can be positive or negative, depending on the values of the numerical coefficients n i . For isolated sources, such corrections scale like the square of the graviton mass m 2 ∼ κM 4 times the "Schwarzschild" radius r s corresponding to the given source, and grow linearly with the distance r. Parametrically, the potential takes the form where φ N is the standard Newtonian potential. Linear theory breaks down at large distances, when the second term is of order unity. It would be interesting to try and match this solution to a non-perturbative exact solution which is well behaved at infinity. As we stated before, the case of no correction to the Newton's law corresponds to the case where (19) is independent of β or γ (cfr. (29)).
Finally, we note that the simple interaction term (2) first considered in [25,28] happens to land on the special case n 2 − n 0 + n 3 = 0, where the above expressions for the gauge invariant potentials are singular. The origin of the singularity is the following. After substitution of the constraints (60), the linearized action no longer depends on ∆E + . In particular, the absence of this variable results in the unwanted restriction (62) on the sources 11 . Nevertheless, beyond the linear order, the action will depend on ∆E + , and hence the "restriction" will no longer exist. Rather, a nonlinear equation will determine the value of ∆E + . Can we nevertheless try to find classical solutions in a perturbative expansion?
The above considerations suggest an expansion scheme for the singular case n 2 − n 0 + n 3 = 0, where E + is treated as a much bigger quantity than the rest of the linearized fields 12 (such as ψ). Taking For distances shorter than the inverse graviton mass, we have ∆E ∼ φ 1/2 N , and hence we may This accidental symmetry is similar to that which exists in ordinary massive gravity where the linear action propagates 5 dof whereas a new ghost-like dof appears at the non-linear level [35,26]. However, in that case the accidental symmetry corresponds to a symmetry of the massless theory and no further constraints are needed in the sources. 12 Some of the linearized fields will be of the order of E as is clear from (60). 13 All the coefficients will have corrections of order O(h). However, for the rest of coefficients one expects that they will yield second order small corrections.
At distances which are large compared with the inverse graviton mass, ∆E ∼ (∆φ N /m 2 ) 1/2 , and we expect These very crude arguments seem to indicate that, also in this special case, there is no vDVZ discontinuity. However, for finite m, there are significant modifications to the value of the "gauge invariant" potential Φ which determines the motion of slowly moving particles. For isolated sources, such modifications scale like r 1/2 s , where r s is the "Schwarzschild" radius corresponding to the given source. They grow with the distance as r 3/2 below the graviton Compton wavelength m −1 , and as r 1/2 for larger distances. The potential Φ becomes of order one for r m −2 r −1 s , beyond which we enter a non-perturbative regime. It would be interesting to confirm this heuristic analysis in a numerical study of a spherically symmetric solution with sources. This is left for further research.

Perturbation theory of Proportional de Sitter Metrics
As stated in Section 2, another interesting class of solutions of bigravity can be constructed from two proportional metrics with a constant proportionality factor. Let us define our perturbations as All indices will be handled with the Ω µν metric.
We first focus in the interaction term for a general potential (4). Using (21) we can writẽ where indices are manipulated with the metric Ω µν , e.g. h g = Ω µν h g µν , and We have also introduced an effective Newtons's constant κ + for later convenience. Note that the massive graviton corresponds to h + µν = (h g + h f ) µν . This is to be expected, as for h g µν = −h f µν the metrics are still proportional and therefore the perturbations are standard massless gravitons of GR in vacuum. Also, in the present set-up, h + µν are the quantities invariant under the diagonal diffeomorphisms. Notice also that the mass term does not have in general a Pauli-Fierz form, This particular form can only be achieved by properly tuning the parameters. This is in contrast with other ways of getting massive gravitons, such as dimensional reduction, where the original symmetry group is much larger. Here, the degrees of freedom of the original theory are 8 which can be split into a massless graviton with 2 polarizations and a massive graviton with 6 polarizations 14 .
The expression of the massless graviton as a linear combination of the metric perturbations will be given below. From Eq. (71) we note that whenever m t = 0 there is an enhancement of the gauge symmetry, which now admits all transformations which leave the traces h g and h f invariant 15 . This corresponds to the transverse subgroup of the diffeomorphisms, which has been recently considered in [37]. In this special case the gauge symmetry is enough to have just two massless gravitons propagating 16 .
Let us now consider the case of generic m s and m t . For simplicity we will concentrate on perturbations around de Sitter solutions which will be foliated by spatially flat sections, where a(η) = −(Hη) −1 , H 2 = Λ g /3 being a constant and η ∈ (−∞, 0). The kinetic term in (1) will be given by (cfr. (71)) with Λ f = γ −1 Λ g . To second order in perturbations we can rewrite the kinetic term in terms of a massive and a massless field, where κ − = κg 1+κ , κ + = κ g κ −1 (1 + κ), with κ = γκ g κ −1 f , g −µν = Ω µν + h − µν and g +µν = Ω µν + h + µν . Besides, we have introduced the massive and massless combinations The dynamics of the massless part is well known. One easily finds that only the tensor modes are dynamical. For the generic massive theory in de Sitter space, studying the longitudinal mode of the massive representation we would argue that the only ghost-free possibility is the Fierz-Pauli mass term, m 2 t = −m 2 s [2,41]. However, in general, this mode decouples only at high energies (larger than a combination of the rest of relevant mass scales). For intermediate energy scales, the longitudinal mode is coupled to another scalar mode which can modify this picture [16,27]. Also, the curvature scale H could play a role in making these intermediate scales phenomenologically relevant. We will study this possibility directly in the unitary gauge.
Let us first split the degrees of freedom of the massive combination into scalar, vectorial and tensorial modes, where ψ, B, and E are the scalar modes, F i and V i are vector modes, and t ij is a tensor mode. The vector modes are divergenceless and the tensor modes are transverse and traceless. The expansion of the kinetic term in this foliation can be extracted from the usual expansion in de Sitter space (see e.g. [42], notice however the difference of convention). One finds where H = a(η) ′ /a(η) = a(η)H and the prime refers to derivative with respect to the conformal time η. We have also introduced the d'Alembertian = η µν ∂ µ ∂ ν and the Laplacian △ = ∂ i ∂ i . The interaction term (71) reads We can now analyse the different components in turn.

Tensor and Vector Modes
The action for the massive tensorial modes is simply From this equation we can read the mass of the graviton which will be given by m 2 t , and the tachyon-free condition will simply read m 2 t ≥ 0.
Regarding the vector modes, their action is The field V m enters the action without time derivatives, and thus its variation yields the constraint, Taking this constraint into account, the action for the vectorial modes up to second order can be written as This Lagrangian has the usual signs, and thus no ghost or tachyons appear in the theory for m 2 t ≥ 0. More concretely, we can canonically normalize the previous field equation with the field redefinition We conclude that the only constraint we get form the analysis of the vector and tensor modes is

Scalar Modes
From (79) and (80), the second order Lagrangian for the scalar part reads (2) B is non-dynamical, and for m 2 t = 0 it is determined in terms of the other fields. For m 2 t = m 2 s , A appears only linearly in the mass term. For the flat case H = 0 and a(η) = 1, this makes A a Lagrange multiplier and thus its variation gives rise to a constraint between the fields E and ψ, leaving just one scalar propagating degree of freedom. In the de Sitter case, the result is the same, although this is not so obvious from the previous expression for the action until one substitutes the constraints.
The variation with respect to A and B yields the constraints so the condition is considerably relaxed at wavelengths shorter than the inverse graviton mass.
Once we have established the positivity of part of the Hamiltonian, let us see what happens to rest of it, namely to the potential part. This part will be given by where the coefficients are rather cumbersome and we omit them. Before proceeding, it should be noted that the Hamiltonian we are considering is time dependent, and hence not conserved. Its positivity and boundedness is a useful criterion only as long as we consider time-scales shorter than the expansion time, or energies larger than H. This is what we may call the adiabatic limit.
Hence, let us assume that m s , m t ≫ H, even if their difference is much smaller m 2 s − m 2 t H 2 , so that we can satisfy the positivity of the kinetic term as discussed above. We have checked that within this adiabatic limit, the potential V grows negative and unbounded below for −∆/a 2 ≫ m 2 . Instabilities at high momenta have been previously studied in [43], and they are just as bad as ghost instabilities. Unlike the case of tachyons, the phase space for instability is infinite and this yields infinite decay rates.
If the masses m t and m s are small, of order of the expansion rate H, then we are outside of the adiabatic limit, and the Hamiltonian above is not a very useful indicator of stability. Instead, we should use a conserved charge associated to the time-like Killing vector for length scales smaller than the horizon [44]. Due to the existence of the cosmological scale, it is in principle possible (although by no means clear) that there may be some range (with α > 2), where this conserved charge is positive definite. The effective theory would then be well defined for momenta larger than H (corresponding to modes within the horizon), provided that the theory is cut-off at the energy scale m(H/m) α/2 . We leave the study of this conserved charge for further research. We note, however, that we need a theory which is applicable to wavelenghs much smaller than the horizon −∆/a 2 ≫ H 2 , where the adiabatic approximation should again be valid. We have checked that for −∆/a 2 ≫ H 2 m 2 , the potential V grows negative and unbounded below, so the possibility of a range of the form (99) where the conserved charge is positive does not look particularly promising.
Finally, for the case m 2 s = m 2 t the analysis of the degrees of freedom has already been performed in another foliation in [45] (see also [46]). In our analysis for this case we find M 2 (η) = M 3 (η) = 0 and thus ψ is not a propagating field. After varying the action with respect to ψ we obtain a constraint which after substitution yields the Lagrangian This Lagrangian will be ghost-free and tachyon-free for µ H ≤ 0. This reduces to the well known condition m 2 ≥ 2H 2 [47]. Another case which differs from the usual Fierz-Pauli term and is still well defined is the case of Lorentz-breaking mass terms [15,16]. These mass terms arise from solutions of the form sketched in Section 2. A simple example would be given by two de Sitter metrics with λ i = {βγ −1 , γ −1 , γ −1 , γ −1 }. The interesting ranges propagating only one scalar or no scalars at all are the same as those in [15,16], and for generic mass terms with two scalars, there is always a gradient instability at high energies as happens in the Lorentz preserving case [39].

Conclusions
Bigravity is a natural arena where a non-linear theory of massive gravity can be formulated. This theory is invariant under the "diagonal" group of diffeomorphisms, under which both metrics transform.
We have shown that if the interaction between the two metrics is non-derivative, we can always find exact solutions in the Schwarzschild-(A)dS family, with Lorentz breaking asymptotics. These are "Type I" solutions [28], in which both metrics cannot (generically) be brought to diagonal form in the same coordinate system. For a given interaction potential, the degree of Lorentz-breaking is an adjustable parameter.
Perturbations around Lorentz-breaking bi-flat (or bi-de Sitter [39]) solutions lead to gravitons with Lorentz-breaking mass terms. Because of the invariance under diagonal diffeomorphisms [22], mass terms with components h 0i are absent from the second order Lagrangian. This, in turn, leads to a well behaved theory of linearized perturbations [22,31], which is not afflicted by the vDVZ discontinuity. It is somewhat puzzling that in the linear theory, there are corrections to the Newtonian potential which are proportional to the square of the graviton mass and which grow linearly with the distance to the origin. On the other hand, as mentioned above, these theories admit the Schwarzschild metric as an exact solution for the same values of the parameters. Thus, the linearized solutions for a static spherically symmetric sources do not coincide with the linearization of the known vacuum solutions. This seems to indicate that this theory has a linearization instability such as the one which is found in other contexts [48], some of which are related to massive gravity and may have important phenomenological consequences [49]. Another possibility is that there may be other exact solutions which coincide with the linearized approximation at large distances, and those may be the relevant ones which can be matched to spherically symmetric matter sources near the origin. This issue clearly deserves further investigation.
We have also considered perturbations to solutions where both metrics are proportional to each other, focusing in the case of de Sitter. This has led us to consider generic mass terms for gravitons in de Sitter space, beyond the Fierz-Pauli case. We find that for the case of Lorentzinvariant mass terms, only the Fierz-Pauli combination is free from instabilities at high momenta −∆/a 2 ≫ m 2 (and only for m 2 ≥ 2H 2 ). For the Lorentz-breaking case, the situation is analogous to that in flat space [15,16,39].