Heavy spin-2 Dark Matter

We provide further details on a recent proposal addressing the nature of the dark sectors in cosmology and demonstrate that all current observations related to Dark Matter can be explained by the presence of a heavy spin-2 particle. Massive spin-2 fields and their gravitational interactions are uniquely described by ghost-free bimetric theory, which is a minimal and natural extension of General Relativity. In this setup, the largeness of the physical Planck mass is naturally related to extremely weak couplings of the heavy spin-2 field to baryonic matter and therefore explains the absence of signals in experiments dedicated to Dark Matter searches. It also ensures the phenomenological viability of our model as we confirm by comparing it with cosmological and local tests of gravity. At the same time, the spin-2 field possesses standard gravitational interactions and it decays universally into all Standard Model fields but not into massless gravitons. Matching the measured DM abundance together with the requirement of stability constrains the spin-2 mass to be in the 1 to 100 TeV range.


Introduction
Numerous cosmological and astrophysical observations have confirmed the presence of a Dark Matter (DM) component in our Universe. Until now this unknown type of matter has been seen only through its gravitational interactions, which resemble those of ordinary matter. Its effects are visible in the rotation curves and velocity dispersions of galaxies, gravitational lensing, matter distribution power spectra, structure formation, Baryon Acoustic Oscillations and angular power spectrum of the Cosmic Microwave Background [1]. The standard paradigm treats the unknown DM particle as a cold relic density which has been created through a model-dependent production mechanism in the early Universe. General Relativity (GR) as the theory for gravity (including a cosmological constant Λ which accounts for the observed amount of Dark Energy) together with a particle physics model for cold Dark Matter (CDM) yield the concordance description of cosmology, the ΛCDM model. See [2] for a recent review of its status quo. The most popular DM models moreover assume that the DM particle is weakly coupled to baryonic matter and hence might be produced in colliders, directly detected in dedicated experiments or indirectly observed through astro-particle signatures. From a theoretical perspective, many of these models lose some of their attractiveness because, typically, they are either not very well motivated from fundamental principles or they introduce a large number of unobserved additional fields (such as Supersymmetry). Unfortunately, on the experimental side, all attempts to produce or detect the DM particle have remained unsuccessful so far [1,[3][4][5][6][7].
The absence of any signatures for DM apart from its gravitational effects motivates a shift of paradigm in the way to think about the nature of DM. Instead of augmenting the Standard Model (SM) by an additional field, we suggest that the DM particle may instead arise in a minimal extension of the gravitational sector, namely in the form of an additional massive spin-2 field. To us this seems to be a natural and well-motivated proposal, since there is no evidence supporting the fact that DM shares the quantum numbers of one of the SM particles and we only observe it through its gravitational interactions.
General Relativity can be treated as the unique theory of a single massless spin-2 particle, the graviton. We will colloquially refer to this point of view as the standard description of gravity. Since massless spin-2 fields cannot interact with each other [8], the most natural and minimal addition to gravity is that of a massive spin-2 field. Studying the effects of a massive spin-2 field in addition to standard gravity amounts to answering fundamental questions of field theory. Remarkably, for several decades it was believed that no consistent theory for gravitating massive spin-2 fields can be formulated owing to the unavoidable presence of a fatal ghost instability [9]. Only recently has the unique description which avoids the ghost been found [10][11][12][13]. Since it involves an additional dynamical tensor field, which mixes with the gravitational metric, the corresponding theory has been titled "bimetric theory of gravity". If fundamental massive spin-2 particles exist, they are described by this unique theory, which automatically leads to a modification of gravity. For the history and detailed reviews of theories for massive spin-2 fields we refer the reader to [14][15][16].
Following ideas outlined in [16,17], it has been proposed that the existence of a massive spin-2 particle can explain all the effects related to DM [18] (see also [19]). The present work is dedicated to providing details of and more insights into this novel proposal.
Summary of results. Being a modification of gravity, bimetric theory must satisfy constraints coming from Solar System tests and cosmology. We confirm that a large value for the spin-2 mass, together with a small value for the "second Planck mass" of the metric that does not couple directly to matter, imply that the static spherically symmetric and cosmological solutions of bimetric theory always resemble those of GR. In this parameter region, where bimetric theory passes all observational tests of GR, the additional massive spin-2 field continues to gravitate but decouples from matter, automatically providing an ideal DM candidate.
We derive the conditions which ensure the validity of a perturbative treatment of bimetric theory for the interesting parameter regions and energy regimes. The structure of cubic and higher interactions for the spin-2 fields forbids a decay of the massive field into massless gravitons, resulting in a discriminating feature of the bimetric model.
Requiring sufficient spin-2 DM to be produced in the early Universe and imposing constraints coming from its possible decay into SM fields, we obtain the allowed region in the bimetric parameter space: The spin-2 mass has to lie within the narrow region of 1 TeV m FP 66 TeV and the ratio of Planck masses must satisfy 10 −11 α 10 −15 . This region overlaps with the one where classical solutions to the bimetric equations resemble GR to a very high precision and therefore the theory passes all observational tests. Moreover, our setup introduces no additional energy scale significantly higher than the weak scale and thus does not create any new hierarchy problems with respect to GR.
Our novel DM proposal naturally explains the absence of signals in (in)direct detection experiments and colliders. A prediction of the model is that any future experiments of this kind will continue to produce null-results. In turn, a detection of a DM particle with mass below our predicted value would rule out our proposal of a heavy spin-2 as sole explanation for the observed DM. Alternative tests of our proposal could be based on its gravitational or its self-interacting nature.
Outline of the paper. Section 2 is dedicated to reviewing relevant details of bimetric theory. We provide its action, equations of motion and present the maximally symmetric solutions with the corresponding mass spectrum. In section 3, we discuss the two parameter regions for which the classical solutions of the theory resemble those of GR. The regime of parameters and energies where the bimetric action can be treated perturbatively is derived in the beginning of section 4. Thereafter, we compute the cubic and quartic vertices in the spin-2 sector, verifying the absence of decay terms into massless gravitons. The phenomenology of spin-2 DM is explored in section 5. Finally, we discuss our results in section 6. Additional supporting details can be found in the appendices.

Action and equations of motion
In order to set down some notation and facilitate our later discussions we first provide some of the required basic details of the ghost-free bimetric theory. For further details and a recent review on the subject we refer to [16]. The theory is defined by the action [13], |g|R(g) + α 2 |f |R(f ) − 2α 2 m 2 g |g| V (S; β n ) + d 4 x |g| L m (g, Φ) . (2.1) Here m g is a mass scale related to the reduced Planck mass via, The dimensionless α 2 (the "ratio of Planck masses") measures the relative interaction strength of the two tensor fields. 1 In addition, the interaction potential V (S; β n ) contains 5 dimensionless parameters β n . Of these, β 0 and β 4 act as bare cosmological constants for g µν and f µν , respectively, and therefore encode nonlinear self-interactions while the remaining parameters β 1 , β 2 , β 3 encode the nonlinear interactions between the two tensor fields. 2 The form of the interaction potential V is constrained by demanding absence of the so called Boulware-Deser ghost [9] and is given by [10,13], where the e n (S) are the elementary symmetric polynomials defined in terms of the eigenvalues of the matrix S. Explicitly they can be obtained via tracing the unit weight totally antisymmetric products, e n (S) = S µ 1 [µ 1 · · · S µn µn] . (2.4) From the definitions it follows that e 4 (S) = det(S) and that e n (S) = 0 for all n > 4. The matrix argument S which appears in the interaction potential is a square-root matrix defined through the relation, For invertible S, the identity e n (S −1 ) = e 4−n (S)/e 4 (S) can be used to show that, 1 As we will see later, α also quantifies the mixing between the propagating states relative to the interacting states. In the literature it is common to use mass parameters mg and m f in front of the kinetic terms as well as a mass parameter m in front of the interaction terms. The scale of m is redundant since it only parameterises the overall scale of the βn and can be chosen freely without loss of generality. In our notation, the relation to these additional mass parameters is given by m f = α mg and m = √ mgm f .
Apart from the matter interactions, the structure of the theory is therefore completely symmetric in terms of g µν and f µν . In fact, in the absence of matter couplings the structure of the bimetric action is invariant under the following discrete interchanges, This property is quite useful since it allows us to obtain the f µν equations of motion directly from the g µν equations.
Finally, the matter Lagrangian L m contains the SM matter fields Φ, which we have taken to be minimally coupled only to g µν here. This choice is without loss of generality since the theory treats the metrics symmetrically and, in the bimetric theory, standard matter fields can only couple minimally to one of the tensor fields without introducing ghost instabilities [20,21].
The equations of motion that follow from the action (2.1) are given by, Here G µν (g) = R µν (g) − 1 2 g µν R(g) is the Einstein tensor computed with respect to g µν and is the Einstein tensor computed with respect to f µν . The interaction contributions V µν ,Ṽ µν and the matter stress-energy T µν are defined by, As noted, (2.8b) can be obtained directly from (2.8a) by making use of the interchange symmetry (2.7). The interaction contributions are matrix polynomials in S, which can be written, where the tensors Y (n) (S) are defined as, For example, written out explicitly we have that, (2.12) We note that V µν andṼ µν as written in (2.10) are symmetric in their indices, although not manifestly so. This follows from the fact that both S and S −1 are symmetric whenever their indices are raised or lowered using either of g µν or f µν . 3 The theory defined by the action (2.1) is generally covariant under the diagonal group of common diffeomorphisms. The fact that the interaction potential is covariant on its own implies the following divergence identities (see e.g. [24]), as well as the algebraic identities [25] (see also [26,27]), where V is the interaction potential (2.3) appearing in the action (2.1). For a covariantly conserved source, the standard Bianchi identities, ∇ µ G µν = 0 and∇ µ G µν = 0, imply the constraint equations ∇ µ V µν = 0 and∇ µṼ µν = 0. Due to the identity (2.13), these are not independent and so in all only give 4 constraints. Apart from that, an additional scalar constraint can be constructed [28] (which was first found in the Hamiltonian formulation [13]). These 4 + 1 constraints serve to remove 5 dynamical modes from the 10 + 10 = 20 components of the two tensor fields. The diffeomorphism invariance removes 2 × 4 = 8 more. Bimetric theory therefore propagates 20 − 8 − 4 − 1 = 7 degrees of freedom. As we will see next, when such a split makes physical sense, these degrees of freedom correspond to a massless spin-2 field (2) and a massive spin-2 field (5).

Proportional solutions & mass spectrum
An important class of solutions in bimetric theory without any matter sources are the proportional solutions, defined byf µν = c 2ḡ µν . For such an ansatz the Bianchi constraints immediately imply that c 2 is a constant. In order to simplify notation we will set c 2 = 1 in what follows. This can be done without any loss of generality by scaling f µν and properly redefining α along with the β n . Such a scaling is possible since we do not couple f µν to matter in our considerations, which gives rise to a redundancy in the parameter space.
For the proportional ansatz the bimetric vacuum equations reduce to [22], with constants, Consistency between the equations now requires Λ g = Λ f ≡ Λ. Since we have set c = 1 this relation generically fixes one of the β n parameters. 4 This class of solutions thus corresponds to the maximally symmetric solutions of GR. Flat space solutions with Λ = 0 require fixing one of the β n and whenever we discuss flat backgrounds this will always be implicitly assumed. For spacetimes admitting Poincaré or (Anti) de Sitter isometries the representation theory of spin-2 fields is well known. It is therefore natural to study perturbations of the proportional solutions. Perturbation theory in bimetric theory is notoriously challenging due to the presence of the square root matrix in the interaction potential and the general problem was only recently resolved [28][29][30]. For the proportional solutions, however, the situation simplifies greatly.
We define linear fluctuations h, around the proportional backgrounds by, The canonically normalised mass eigenstates are then defined through [22], where, for future reference, we also note the inverse relations, The parameter α thus quantifies the mixing between the fluctuations. In terms of the mass eigenstates (2.18) the quadratic part of the action (2.1) diagonalises into (indices are raised and lowered usingḡ µν ), where L GR is the quadratic theory obtained from the Einstein-Hilbert action including a cosmological constant, i.e. |g|(R−2Λ). The detailed expression for this is given in eq. (B.7). We have defined the Fierz-Pauli mass of the massive spin-2 field, (2.21) Note that our parametrisation implies that the parameters β 1 , β 2 , β 3 are on the order of m 2 FP /m 2 Pl . The quadratic theory contains a massless graviton δG µν , which mediates standard gravitational interactions with Planck mass m Pl and an additional massive spin-2 field δM µν with mass m FP . Note that the massive spin-2 field couples to the matter stress-energy and therefore also mediates gravitational interactions but with a coupling α/m Pl .
For small α, the matter coupling of the massive field will be suppressed with respect to that of the massless field. One may therefore expect to recover a situation close to linearised GR for small enough α. On the other hand, as we will see later when studying higher-order interactions, while the massive mode decouples from the SM matter, it does not decouple from gravity in the α → 0 limit. In fact, it continues to gravitate with the exact same strength as SM matter. This makes the massive spin-2 field of bimetric theory an interesting candidate for a DM particle. Similarly, another way to recover linearised GR (at least at low energies) is to consider large values for m FP , which also decouples the heavy spin-2 field from the matter sector. The following section is dedicated to a detailed discussion of the behaviour of the theory in these two parameter regimes.

Recovering General Relativity
As any other modification of GR, bimetric theory containing a second tensor field generically changes the laws of gravity. Since GR is well-tested over a large range of energy regimes, we need to carefully evaluate the observational constraints on bimetric theory and make sure that its predictions do not differ too much from GR. In this section we will see that there are two different (but overlapping) regions in the parameter space of bimetric theory for which certain classical solutions for the physical metric approach those of GR. In particular, the cosmological as well as the static spherically symmetric solutions to the bimetric equations of motion both resemble GR in the overlap of these two regions.

The GR regimes for physical solutions
The two separate parameter regions which recover GR for the physical metric g µν can be motivated based on the linear theory around proportional backgrounds: (i) The more general option is to consider a large hierarchy between the "Planck masses" of the two metrics, i.e. α 1. Physically this corresponds to a very feeble coupling of the massive spin-2 field to matter sources, irrespective of its mass. It also implies enhanced self-interactions of the massive field and a large value for the physical Planck mass. All known solutions of bimetric theory coincide with GR solutions for g µν in the limit α → 0.
(ii) The second option is to take the Fierz-Pauli mass m FP to be large, typically m 2 FP Λ, which effectively means that ξ in (2.21) should obey Λ/m 2 Pl ξ 2 . Additionally, since ξ sets the scale of m FP in units of m Pl , we should also require ξ 1. Regarding the massive spin-2 field as DM will turn out to give much more stringent bounds. In physical terms, we would like to make the massive spin-2 field heavy enough such that it effectively decouples from the low-energy theory. This option presumably recovers GR for g µν only in the linear regime around the proportional backgrounds, since the notion of m FP has no clear meaning away from these solutions. Nevertheless, this criterion turns out to be useful also in the context of cosmological solutions.
Our analysis in section 5 will reveal that observations favour a combination of both these options. We will therefore discuss the two parameter regions in more detail for two physically important classes of solutions, the static point-source and the cosmological solutions.

Static spherically symmetric solutions
Local gravity tests tell us that any theory for gravity inside the Solar System, up to 10 µm, must follow the predictions of GR to high precision [31] and bimetric theory studied in the context of DM has to pass these tests. To approximate modified gravity effects inside the Solar System, one considers static spherically symmetric solutions around a massive source (which would correspond to the Sun). The gravitational field computed in this approximation must effectively resemble GR up to the precision available so far. Another important aspect of studying spherically symmetric solutions is to fix the value of Newton's constant. In modified gravity theories, the value derived in this way can, in principle, differ from the corresponding value obtained in cosmology. Upon comparing local and cosmological observations, this may lead to extra constraints on the theory. Historically, the first attempt to build a massive gravity theory -Fierz-Pauli massive gravity -has been rejected precisely because it fails even basic Solar System tests. The problem arises because the spin-0 mode of the massive graviton adds an extra (fifth) force to gravitational interactions and does not decouple in the limit of small graviton mass, m FP → 0. This effect is known as vDVZ discontinuity [32,33] and we review it briefly in appendix A. It was conjectured in [34] (and confirmed explicitly much later [35][36][37], see also [38]) that the inclusion of nonlinear interactions for the spin-2 field cures this problem and that GR is restored in the limit of small graviton mass. Today this feature is known as the Vainshtein mechanism and it operates in a plethora of modified gravity models, see e.g. the review [39].
As we will see below, bimetric theory with a very heavy spin-2 mass does not require the Vainshtein mechanism, since the solution is linear all the way down to very small lengths, where gravity is not yet tested. In this case, the spherically symmetric solutions recover GR despite being linear, which is in sharp contrast to massive gravity, for which the linear regime always leads to contradictions with Solar System tests. This is a consequence of the fact that massive gravity contains only one propagating massive graviton, while bimetric theory has an additional massless graviton. On the other hand, for bimetric parameters which require us to go beyond the linearised approximation, we still have to rely on the Vainshtein mechanism to restore GR. In this case, the restoration through the Vainshtein mechanism is quite similar to massive gravity.

Derivation
An appropriate ansatz for spherically symmetric solutions in bimetric theory reads [40], ds 2 (g) = −e ν dt 2 + e λ dr 2 + r 2 dΩ 2 , (3.1) where ν, λ,ν,λ and µ are functions of r. We will always assume that the functions (ν, λ,ν,λ) are much smaller than unity, corresponding to weak sources. On the other hand, the function µ can be either small or large and in the latter case carries information about nonlinear effects. In massive gravity the function µ can be associated with a Stückelberg field [37]. The common diffeomorphism invariance has been used in the above ansatz to remove a function in front of r 2 dΩ 2 in the g µν metric.
Linearising the equations of motion (2.8), one can obtain the following solutions [40] (see also [41,42]), where C 1 and C 2 are two integration constants, to be fixed by matching the solution to the source. Depending on the parameters of the model and the mass of the central source, the linearised approximation may not be valid for all distances. The above expressions were obtained under the assumption that nonlinearities in µ can be neglected. The Vainshtein mechanism starts to operate exactly when nonlinearities in µ become important. The equations of motion (2.8) can be solved analytically in a different regime, which does not rely on linearity in µ. Assuming that we are deep inside the Compton wavelength, r m −1 FP , one finds a seventh order algebraic equation for µ [40]. The solution in this regime is valid down to small radii, and it can be matched to a solution inside the source. To this end, we introduce the Schwarzschild radius, where ρ is the density inside the central source and R is the radius of the body. Note that the above expression has an extra factor (1 + α 2 ), resulting in an extra factor (1 + α 2 ) in Newton's constant with respect to GR. Another relevant scale is the Vainshtein radius, below which the nonlinearities in µ kick in. One then finds that for r V At smaller radii, r r V the solution changes its form to, which restores GR for the physical metric g µν . This is precisely the Vainshtein mechanism operating for radii r r V . The constant expression for µ depends on the parameters of the Lagrangian as well as the exact form ofλ andν.
The matching of the linearised solution (3.3) and the solution inside the Compton wavelength (3.6) fixes the constants of integration C 1 and C 2 as follows, Note that in the linear regime, the metric functions of g µν receive an extra factor 1/(1 + α 2 ) with respect to their behaviour in the Vainshtein regime. Moreover, the above matching is only valid when the Vainshtein regime is present at all. Otherwise, when the linear regime is valid all the way down to the source, the solutions obtained in the linearised approximation must be matched to the source. We discuss this case below, in the context of large values for the spin-2 mass. It is worth pointing out that the scale of nonlinearity r V is not directly related to the validity of the perturbative expansion for bimetric theory which we will address in section 4.1.2. The scale where classical solutions become nonlinear depends on an extra scale of the problem, namely the mass of the central source.

Λ
Let us now discuss the region where the Fierz-Pauli mass is large, i.e. m 2 FP Λ. In the limit of infinitely large mass, both the Compton wavelength m −1 FP and the Vainshtein radius r V vanish. This means that the linearised approximation is valid for all radii. From (3.3) one can easily see that for m FP → ∞ we find the GR solution, Note that, in this case, one cannot use the expressions in (3.8) for the integration constants since they were obtained by assuming that the Vainshtein regime operates for small radii. In the limit m FP → ∞ the solution always remains linear and one needs to redo the matching to the source. We will not go into the details of this computation, but only give the result.
Assuming that λ andλ are given by (3.9) outside the source, and by the same expressions, but with C 1 being a function of the radius inside the source, one obtains, (3.10) We conclude that the local Planck mass coincides with our original definition of m Pl . This result could have been anticipated from the action written in the terms of mass eigenstates (2.20). Since the linear approximation is valid in the limit m FP → ∞ (at least outside the sources), the quadratic action (2.20) is sufficient for studying spherically symmetric solutions. Note that both the massless δG and massive δM spin-2 field contribute to the physical metric g µν via the relation (2.19a). The coupling constant between the massless graviton δG µν and the source is precisely m −1 Pl . At the same time, the massive spin-2 mode δM µν is also excited by the source term, but because of the vanishing Compton wavelength m −1 FP , the solution for δM µν outside the source vanishes in the limit m FP → ∞. As a consequence, the only contribution to the physical metric h µν comes from the massless mode δG µν , c.f. (2.19a). Hence we recover exactly the results in (3.9) and (3.10).

The region α 1
Now we turn to the limit α → 0, where at the same time we keep m FP constant. Neither the Vainshtein radius in (3.5) nor the Compton wavelength m −1 FP vanish in this limit. Therefore, in contrast to the case m FP → ∞ considered above, the theory enters a nonlinear regime for small enough distances (at least for large enough r S ). According to the general discussion in section 3.2.1, the solution in the linear regime is valid for r r V . From (3.3) in the α → 0 limit we then find, where we also used (3.8). This shows that, for large radii, GR is restored. The physical Planck mass is again m Pl since from (3.4) with α → 0 we get, Notice that, in contrast to the case of large mass, the function µ (the "Stückelberg field") does not vanish, At r ∼ r V it becomes of order unity, confirming that the linear approximation breaks down at this scale. Thus we have to resort to the nonlinear regime for r < r V . Using (3.6) and (3.7), it is straightforward to show that, in the limit α → 0, the solution in the nonlinear regime with r < r V is again given by (3.11) and thus coincides with GR.
Once more we could have started from the quadratic action (2.20) and anticipated part of the result. Solving the equations for the massless eigenstate δG µν , we recover the GR solution as long as the linear regime is valid. Newton's constant is given in terms of the Planck mass m Pl , as can be read off from (2.20). The relation (2.19a) shows that the only contribution to h µν is δG µν in the limit α → 0 and hence h µν in the linear regime corresponds to a GR solution. This is in accordance with (3.11), which was obtained from the general formalism of spherically symmetric solutions. On the other hand, when nonlinearities kick in, we cannot rely anymore on the quadratic action (2.20) since the nonlinear terms become as important as the linear ones. In other words, (2.20) is not sufficient for studying the behaviour of the metric inside the Vainshtein radius. In this regime GR is restored by nonlinear effects, c.f. (3.7), which is independent of taking any parameter limit.  Figure 1. Schematic diagram showing how GR is restored for different distance scales r, depending on the spin-2 mass m FP . In the red-shaded region, the solution is nonlinear and GR is recovered via the Vainshtein mechanism. In the linear regime (blue-striped region), for r < m −1 FP , it is necessary to require α 1 in order to recover GR (white region), while for r m −1 FP , GR is recovered due to the exponential fall-off of the Yukawa potential (green-shaded region). We conclude that in the limit of small α, GR is restored for all radii, with Newton's constant given by the Planck mass m Pl . It is worth emphasising that this type of GR restoration for the physical metric g µν is quite nontrivial, since it involves a transition between the linear and nonlinear regimes. All features defining the regime of the solution are hidden in the second metric f µν , including the Stückelberg field µ. The different regimes in which GR is recovered for static spherically symmetric solutions are visualised in Fig. 1.

Cosmological solutions
Just as the local gravity tests, the GR based ΛCDM concordance model has been confirmed to high precision and therefore puts stringent constraint on modifications of gravity. The homogeneous and isotropic solutions to the bimetric equations of motion were first derived in Ref. [43][44][45]. For general parameters they give rise cosmological observables which differ significantly from GR predictions. The behaviour of cosmological solutions for large m FP has previously been discussed in [46] and for small α in [47]. Here we will take a slightly different approach with respect to these references, which will allow us to treat both cases simultaneously and to show that bimetric theory again resembles GR in the overlap of these parameter regions. For a review of bimetric cosmology, we refer the reader to [48].

Derivation
Restricting our analysis to the bidiagonal case, we can put the metrics on the form, where k = 0, ±1 corresponds to a flat, open and closed universe respectively. It follows from the equations of motion that we must have the same k in both metrics. Here all isometries have been used to put g µν on the standard FLRW form. The "lapse" function X(t) can be solved for directly from the Bianchi constraint to give X =Ẏ /ȧ. The matter source coupled to g µν is taken to be a perfect fluid, T µ ν = diag(−ρ, p, p, p). The spatial scale factors a(t) of g µν and Y (t) of f µν are then solutions to the dynamical equation, with r(t) ≡ Y /a, 16) and the algebraic constraint, The first of these is a modified Friedmann equation for the physical scale factor a(t) while the second determines r(t) in terms of ρ(t). In addition to the linear dependence on the energy density ρ, the squared Hubble function H 2 = (ȧ/a) 2 is now also sourced by the contribution coming from the interaction potential. In general, one solves the polynomial equation (3.17) for r in terms of ρ and plugs the solution back into (3.16). This results in an equation of the form , where F is an analytic function of the energy density whose precise form is determined by the choice of bimetric parameters. Generically, the nonlinear nature of F [ρ] leads to significant deviations of bimetric cosmology from ΛCDM.
In addition to the above equations, the matter source is subject to the same continuity equation as in GR. The various source components in ρ(t) therefore dilute in time in the standard way. As matter dilutes and ρ(t) → 0, we see from the algebraic equation (3.17) that r(t) → const. and thus the metrics become proportional with r 2 being the constant of proportionality. Moreover, it follows from (3.16) that this late-time de Sitter attractor solution has a cosmological constant given by, Here we have taken the asymptotic constant value of r to be r = 1 without any loss of generality (c.f. our discussion in beginning of section 2.2). The fact that there is a de Sitter attractor at late times is crucial for our logic of extracting phenomenology out of these solutions since observations suggest that we are now living in a cosmological epoch dominated by a cosmological constant. From the algebraic equation (3.17) we can immediately infer that as soon as the matter density obeys ρ(t) β i m 4 Pl we can neglect its contribution and solve for r = const. up to small corrections. To quantify these corrections we make a power series ansatz, 5 which we plug into (3.17) and subsequently solve for the coefficients a n . After plugging the resulting expression for r(ρ) back into (3.16), the modified Friedmann equation reads, 6 Before discussing the validity of this expansion and what it implies for the parameters, we make some general remarks. First of all, we note that all but the first correction vanish for β 1 = 0 = β 3 . As can be seen directly from (3.17), in this case it is easy to obtain an exact solution which effectively result in a cosmological constant and a modified Planck mass, matching the first two terms in the right-hand side above. Secondly, we note that the 5 It follows from the analytic implicit function theorem that r(ρ) is indeed analytic around the solution with r = const., provided that 2Λ − 3m 2 FP = 0. 6 In [49] a similar expansion was considered in a bimetric setup with (twin) matter fields coupled also to fµν .
That work only considered the first correction to the Friedmann equation, i.e. the correction to the constant multiplying the term linear in ρ. expansion breaks down when 3m 2 FP = 2Λ. This value saturates the so-called Higuchi bound, 3m 2 FP ≥ 2Λ, which is a well-known unitarity bound for massive spin-2 fields propagating in de Sitter spacetime [50,51]. At the point of saturation a linear gauge symmetry makes the helicity-0 mode of the massive spin-2 field non-dynamical [52,53]. The quest for a nonlinear realisation of this linear gauge symmetry has received a lot of attention lately and the phenomenon has also been studied within the bimetric framework [54][55][56][57]. Although this is a very interesting point in the bimetric parameter space, here we will mainly focus on the regime m 2 FP Λ. As for the validity of the expansion, we note that at nth order in ρ the most dominant correction to the GR term ρ/3m 2 Pl comes with a prefactor on the order of, In order to see this, recall the definition of the Fierz-Pauli mass (2.21) which implies that β n m 2 Pl ∼ m 2 FP for n = 1, 2, 3. Less relevant contributions to the nth order in ρ are suppressed by additional factors of α 2 and/or by powers of Λ/m 2

FP
We conclude that a higher order term in ρ is generically smaller than a lower order term if ρ m 2 Pl m 2 FP . At the present cosmological epoch this is of course quite easy to satisfy even for a tiny mass, since presently ρ/m 2 Λ) we may safely use the expansion to estimate deviations from GR. In order to get a rough order of magnitude estimate for when the validity may break down, we use the fact that at early times we may relate the energy density to the temperature via ρ ≈ T 4 . The bound then implies validity of the expansion for temperatures T 10 9 GeV × (m FP /GeV) 1/2 . In practice, this bound will be even less stringent due to the additional α 2 suppression.

3.3.2
The region m 2 FP Λ Interestingly, despite being physically well-defined only for the proportional solutions, the parameter combinations Λ and m FP play a major role in the expansion (3.20). We see that, apart from the pure cosmological constant term, Λ always enters via the dimensionless ratio Λ/m 2 FP . In particular, the lowest order correction to GR in (3.20) comes as a renormalisation of the physical Planck mass m Pl and is proportional to α 2 Λ/m 2 FP . The strongest bounds on the value of the Planck mass obtained from cosmological/large scale considerations comes from Big Bang Nucleosynthesis (BBN) observations (see e.g. [58]). However, these constrain the value of the physical Planck mass only to within about 10%, or at best a few percent. Hence, this constraint alone does not require a very large m FP or a small α, since the combination α 2 Λ/m 2 FP only has to be less than ∼ 0.1. A stronger motivation for considering m 2 FP Λ comes from the work [46], which considered perturbations of the cosmological solutions discussed above. It turns out that in general a gradient instability is present in the scalar sector which threatens to invalidate linear perturbation theory [59][60][61]. This instability however disappears when (1 + α 2 )Λ m 2 FP and (1+α 2 )ρ 2α 2 m 4 Pl , provided a mild bound on the parameters is satisfied (to wit, β 2 +β 3 > 0). Therefore, considering a large mass is a safe way of ensuring that standard techniques of perturbation theory are still applicable. From the expansion (3.20) it is also clear that the condition m 2 FP Λ alone does not affect all of the corrections to pure GR like behaviour and therefore it cannot serve to fully recover GR from bimetric theory. However, taken together with α 1, deviations of the cosmological solutions from GR are small at all orders in the expansion. As we will see later, if we treat the massive spin-2 field as a DM candidate, then phenomenology indeed favours the region where m 2 FP Λ and α 1. This ensures the compatibility of our model with cosmological observations.
In the literature on bimetric cosmology it is customary to consider a very small mass for the massive spin-2 field, namely m 2 FP ∼ Λ. This is contrary to our approach, but let us briefly comment on the two main motivations for considering a small mass. One is partly historical, relying on intuition from massive gravity where the Vainshtein mechanism is responsible for recovery of GR like behaviour. As we saw in section 3.2, this reasoning is not valid within bimetric theory since, in fact, GR is recovered for large values of m 2 FP without invoking the Vainshtein mechanism. The second motivation is that a small value of the spin-2 mass could lead to a small self-acceleration scale which is "technically natural". This is based on the argument that a vanishing mass restores the full diffeomorphisms and therefore a small value of the mass is protected by a symmetry from receiving large quantum corrections. One may object to this naturalness argument on the ground that i) it perhaps too naïvely carries results from global symmetries over to local symmetries, ii) gravity seems so far to be exempt from obeying any naturalness criterion and iii) there is no fundamental reason to expect naturalness to be a sufficient guide. In addition to these objections, if m 2 FP ∼ Λ, the BBN constraints may actually be more worrisome phenomenologically unless one also requires a small value for α.
We shall not dwell on this issue further but merely note that, if we take the bimetric theory seriously as a model of gravitational interactions, then it certainly seems favourable both from a theoretical and phenomenological perspective to accept a large mass of the additional spin-2 field. Note that this requires us to fine tune the combination of β n parameters in (3.18), in order to produce a small value for the scale of cosmological acceleration.

The region α 1
We stress from the onset that (3.20) is not an expansion in α. Nevertheless, all corrections to the lowest order terms which resemble GR come with at least an α 2 suppression, such that a small α leads to a GR like behaviour of the solution. In other words, although α does not always come into play when comparing higher order terms in the expansion, it does affect the relation between the corrections with respect to the lowest order GR like terms. Note that this is of course already obvious from the nonlinear g µν equations (2.8a) and the modified Friedmann equation (3.16).
In fact, the possibility to restore GR in the quadratic action by taking α → 0 generalises to the full nonlinear level, also beyond the cosmological solutions. To see this, let us recall the bimetric equations of motion given in (2.8). For small α the f µν equations take the form where the O(α 2 ) corrections simply comes from expanding the factor 1/(1+α 2 ) in front ofṼ µν and these can therefore safely be neglected. Now, for regimes where the curvature satisfies R(f ) β i m 2 Pl ∼ m 2 FP , we can neglect also the kinetic term. 7 In this case the f µν equations imply, to first order, that f µν solves the algebraic equationṼ µν =0. The generic solution to the f µν equation for small α is then that the metrics are nearly proportional. As a consequence, the g µν equations assume the form This is consistent with what we found from the cosmological solutions and supports the fact that all solutions for g µν approach GR like solutions for small enough α in the energy regimes where R(f ) β i m 2 Pl ∼ m 2 FP . One may expect that, for small but non-vanishing α, all new effects introduced by the presence of the massive spin-2 mode come in as corrections of O(α 2 ). Interestingly, this turns out not to be the case. For instance, even in the exact α → 0 limit, the bimetric interaction potential contributes to the effective cosmological constant Λ in (3.20) and (3.23), giving rise to background curvature even in the absence of matter. This shows that the universe in bimetric theory can be self-accelerating, i.e. haveä > 0, even in the absence of vacuum energy (i.e. for β 0 = 0).

Heavy spin-2 field coupled to gravity
The arguments of the previous section, which were based on the cosmological and static pointsource solutions, motivate us to further consider the physical implications a heavy spin-2 field coupled to gravity via the ghost-free bimetric interactions. It turns out that such a field naturally has all the desired properties of a suitable DM candidate. In order to elucidate this, we will expand the bimetric action (2.1) in terms of the mass eigenstates defined in (2.18). Before considering the explicit form of this expansion, we discuss some of its general features which are independent of the actual form of the ghost-free interactions. Of course, in the end, we will only consider the specific interactions which are free of the Boulware-Deser ghost.
The discussion of section 4.1 is quite technical and mostly serves to clarify the structure of the expansion and to support its validity. The reader more interested in the final results and their physical interpretation may skip ahead to section 4.2.
7 From the relation (3.18) we infer that some βi may scale as α 2 . In that case the condition on the curvature may turn into R(f ) α 2 m 2 Pl . Later, in section 4, we will indeed restrict ourselves to energies satisfying E α m Pl .

General features
In order to facilitate the subsequent discussion, we first note several useful relations. To this end let us recall the bimetric action (2.1), written here without a matter source and with m 2 = α 2 m 2 g , This gives rise to the vacuum equations of motion (2.8), In what follows we will think of the full nonlinear bimetric action as an infinite expansion around the maximally symmetric background solutions f µν = g µν ≡ḡ µν . It is clear that, absent matter sources, such solutions always exist and imply the background condition Λ g = Λ f = Λ with, The background condition can therefore, in general, be written, The general identity (2.14) can then be seen to imply, So far, this discussion has been completely general and the above expressions hold for any covariant interaction potential V . Let us now restrict our attention to V such that the mass term in the quadratic theory reduces to the Fierz-Pauli one. One can then show that, for fluctuations defined by h µν = g µν −ḡ µν and µν = f µν −ḡ µν , the quadratic theory is always diagonalised in terms of the following canonically normalised mass eigenstates 8 These are of course consistent with our definitions in (2.18a) and (2.18b), which can be seen by using the definition m Pl = m g √ 1 + α 2 .

Mass eigenstates beyond the quadratic expansion
The inverse relations between the metric fluctuations and linear mass eigenstates read, These are valid as linear field redefinitions which diagonalise the quadratic action. When going to higher orders in perturbation theory, one could in principle consider nonlinear corrections to these relations, e.g. add terms of order δG 2 or δM 2 to the right-hand side of (4.7). In other words, beyond quadratic level around the maximally symmetric backgrounds, the definition of mass eigenstates becomes ambiguous. 9 However, in the following we will see that the structure of nonlinear interactions justifies the use of the linear relations (4.7) even at higher orders in perturbation theory. The form of the relations (4.7) implies that we can write the full metrics g µν =ḡ µν + h µν and f µν =ḡ µν + µν as follows, where we have defined a new "background" metric, This structure already hints towards the fact that it may make sense to consider the metric G µν as a massless field which defines the geometry in which δM µν propagates. We now consider the full nonlinear bimetric action as an infinite expansion around the solution g µν = f µν =ḡ µν . We replace the fluctuations h µν and µν around these backgrounds by the linear mass eigenstates using (4.8). Let us start by focussing on the terms involving only G µν , but no δM µν . Neglecting the massive mode, we effectively deal with the nonlinear bimetric action expanded in Clearly, the Einstein-Hilbert terms as well as the β 0 and β 4 terms will just give rise to terms that exactly resemble GR in terms of G µν . The relative factor of α 2 between the kinetic terms serve to give an overall factor m 2 Pl for this part of the action. The only differences, as compared to GR, could come from the potential which is a function of g −1 f . But due to (4.10) we have that the terms without δM µν in this matrix reduce to, and hence (cf. (4.5)), which again is just a cosmological constant term for G µν =ḡ µν + 1 m Pl δG µν . Note that here we have used the background equation Λ g = Λ f = Λ. Thus, the pure nonlinear self-interactions for the linear massless mode δG µν are nothing but the Einstein-Hilbert term, (4.13) In this sense, δG µν behaves exactly like a massless spin-2 field even in its nonlinear selfinteractions.
Next, we turn to the non-minimal couplings to the massive field, where it is particularly interesting to study the terms linear in δM µν . We formally expand the action to linear order in δM µν , treating G µν as a background metric. This is not a valid perturbative expansion in general because we are not keeping all terms of the same order in m −1 Pl . Nevertheless, we can look at the terms linear in δM µν and make formal statements about them which will be valid to all orders in the expansion.
In fact, terms involving δM µν only linearly cancel between the two Einstein-Hilbert terms, √ g R(g) and α 2 √ f R(f ). This directly follows from the fact that the massive fluctuation appears in h µν and µν with a relative factor of α 2 and opposite sign.
For the potential, we simply Taylor expand to linear order in δM in the following way, (4.14) To get to the second line we have used (4.5) and (4.3) together with the definitions of the fluctuations in (4.7). The third line then follows from the background relation Λ g = Λ f . This demonstrates that the bimetric action in vacuum does not contain any terms linear in the massive fluctuations and, in particular, there is no decay into massless gravitons. These general arguments have also been confirmed from explicit calculations of the cubic and quartic interaction vertices using the ghost-free interactions, as we will see in section 4.2. Terms of higher order in δM µν generically do not vanish and give rise to self-interactions of the massive field as well as its non-minimal couplings to the massless graviton δG µν . The two main conclusions of this section are as follows.
(i) The nonlinear self-interactions of the massless eigenstate sum up to the standard Einstein-Hilbert action with a cosmological constant. This is consistent with the interpretation of δG µν as a massless spin-2 field.
(ii) There are no terms linear in δM µν present in the expansion. This implies that there is no decay of δM µν into massless gravitons δG µν at tree level.
These conclusions are very general in the following sense: They are independent of the form of the interactions apart from reproducing Fierz-Pauli theory at the quadratic level and being covariant at the nonlinear level. They are also independent of any nonlinear field redefinitions of the fluctuations. The physical interpretation of the linear mass eigenstates thus seems to make sense even at the nonlinear level and there does not seem to exist any motivation for considering higher-order corrections to (4.7).

Validity of perturbative expansion & absence of strong coupling
We now elaborate in some detail on the validity of the infinite perturbative expansion of the action with α 1 in terms of the mass eigenstates and the related issue of strong coupling. where we have suppressed the index structure along with numerical coefficients. For field values of energy E, these interactions assume the following schematic form, Since here we are only interested in the dependence on E, α and m Pl , we have dropped all the numerical factors. These do not affect our order of magnitude estimates in any dangerous way, but may, together with the tensor structure, at most serve to make some specific combinations vanish. Including them could therefore potentially only serve to sharpen the general and conservative remarks we make here. We also note that, according to the results of the previous section, there will in general be no terms with r + s = 1 present since the bimetric action expanded in mass eigenstates contains no terms linear in δM µν . Nevertheless we will keep these terms for now since they do not influence the general arguments made in this section. For α 1, it is clear that the most suppressed vertices come from pure h µν terms (n = 0) and the most enhanced vertices come from pure µν terms (k = 0). At order m in the fluctuations we thus get vertices with the following structure, where V m denotes a general vertex with m powers of the field fluctuations appearing. Terms with two derivatives, coming from the Einstein-Hilbert terms, will have a similar structure. Their explicit form turns out not to be relevant for the argument of maintaining perturbativity and we therefore focus on the expansion of the interaction potential (but comment on their inclusion towards the end of this section).
It follows directly from (4.17) that in order for this to define a perturbative expansion, we need to require E < α m Pl together with α < 1. In fact, applying the (nth root) Cauchy criterion for convergence show that these requirements are indeed sufficient. Of course, this neglects completely the tensor structure but does provide an estimate for when the expansion is perturbative and the theory is weakly coupled.
A subtlety that arises in ordering the expansion properly is that, for small α, different orders in fluctuations start to mix with each other at some point, depending on their α dependence. For instance, a fourth order vertex may become of the same perturbative order as a cubic vertex which is multiplied by a higher power of α. Therefore, in practice it may be necessary to rearrange the perturbation series differently than in the total number of field fluctuations appearing. This feature is highly dependent on the energy scale under consideration and the exact value of α.
In order to make this more qualitative, let us parameterise the relevant energies as, The condition q > 0 here simply ensures that we only consider energies which satisfy the condition of perturbativity, E < α m Pl for α < 1. With this parameterisation, for a given α, we can probe different energy scales by shifting the value of q. We will return to this point further down in this section. In this parameterisation the above vertices read, In total, the expanded interaction potential therefore consists of vertices V m given by, where V (j) m denotes vertices with m powers of fluctuations and (j + qm) powers of α, with j = 0, 1, . . . , 2m. There is clearly a factor of α suppressing successive terms V (j+1) m with respect to V (j) m . It is however also clear that a term V (j) m+1 , with (m + 1) fields in the vertex, is suppressed with respect to the term V (j) m with m fields by a factor of α q = E αm Pl , which, for q 1, may no longer be small but assume a value rather close to 1. Depending on the exact value of q, the same may be true for V . For any given q 1, we can however always find an integer p ∈ N such that q > 1/p and the number of such terms is therefore always countable.
In order to elucidate this behaviour more closely for small q, i.e. near the perturbativity bound E < α m Pl , let us consider q = 1/p with p ∈ N. It then follows from the definition (4.17) that the vertices V m+jp are always of the exact same order. For definiteness, let us also demand that there is at least a factor of α between dominant terms of what we will refer to as different orders. This means that we consider e.g. V (j) m+k , with k < p, to be regarded of "the same order" as V (j) m . Note that this is of course somewhat arbitrary since the ratio between neighbouring terms in the rearranged series is typically α 1/p . Nevertheless, it allows us to rearrange the perturbation series in the following manner: Here, every line corresponds to terms of "the same order" in the perturbative expansion in the following sense: A line with dominant term V (j) 2 ∼ α j+2/p contains all terms which are of the form α j+(2+k)/p where k = 0, 1, 2, . . . , p − 1 (note however that for a given line the series is not strictly ordered from left to right since, as noted, e.g. V ∼ α j+1+2/p and so on. Therefore lower lines are always more suppressed than the lines above it. It is important to note that every order contains a finite number of terms. Moreover, from the definition in (4.20) we have that V (j) m = 0 for j > 2m. Therefore, for example, the 5 first lines are sufficient for considering the influence of higher order terms on the quadratic vertices and the 7 first lines are sufficient for considering the influence on the quadratic and cubic vertices together and so on.
The kinetic terms have been left out of the discussion so far. They contain two derivatives and give, schematically, where we also have taken into account that the m vertices coming from the Einstein-Hilbert term for f µν will have an additional factor of α 2 in front. It is clear from the structure of (4.22) that the terms can be rearranged in the same way as in (4.21) to give a valid perturbative expansion for E < α m Pl . The kinetic and potential terms will indeed have the same general structure in their expansions. In other words, the structure in (4.21) represents the expansion of the full bimetric action around maximally symmetric backgrounds in terms of the mass eigenstates.
The above discussion illustrates the behaviour of the expansion near the perturbativity bound E < α m Pl . From this we can also understand how the higher order vertices start to influence the physics near these energies. In particular, if we required that none of the higher order vertices play any role at the level of cubic interactions, we would have to restrict ourselves to lower energies. The fact that V (j) m = 0 for j > 2m tell us that the most suppressed cubic term is given by V (6) 3 ∼ α 6+3q . The most dominant quartic term is given by V (0) 4 ∼ α 4q . It is therefore enough to demand that 4q > 6 + 3q, i.e. q > 6, to ensure that all the higher order vertices are subdominant to the quadratic and cubic ones (the cubic are then also automatically subdominant to the quadratic terms). In practice, since we wish to probe energies at least up to the mass scale m FP = ξ m Pl , this means that if we considered ξ = α 1+q with q > 6, we could be sure that no higher order vertices will influence the physics deduced from studying the full cubic theory.
To illustrate this last point with an example, let us consider the minimal mass for which our spin-2 theory can account for all of the observed DM today (c.f. section 5), With this value for α we could be sure that the cubic theory is enough to discuss physics up to energies ∼ 1 TeV and the theory itself remains perturbative up to energies ∼ 0.01 m Pl . As we will see later, our phenomenological analysis reveals that α has to be significantly smaller than 0.01. However, for practical purposes, the perturbativity bound E < α m Pl is sufficient to ensure the validity of all our perturbative expressions in section 5. This is because the inclusion of (a finite number of) higher-order diagrams merely affects the numerical factors of scattering amplitudes (through additional terms multiplied with positive powers of E α m Pl < 1). These corrections are irrelevant for our order-of-magnitude estimates. We will comment on this further below.

Structure of the cubic vertices in ghost-free bimetric theory
The detailed expression for the bimetric action (2.1) expanded up to cubic order in the mass eigenstates is provided in Appendix B. Here we highlight and discuss some of its main characteristics, which also confirms the general arguments given in the previous section.
The dominant structure of the cubic interactions is summarised in Table 1. This table displays the overall coefficients of the various cubic interaction vertices in the δG, δM action. The first column shows the self-interactions of the massless field. As can be verified from the exact expressions in Appendix B.2 these self-interactions are exactly of standard GR form. This is indeed also true for the quartic vertices, i.e. the δG 4 terms exactly match the Einstein-Hilbert structure (c.f. Table 2). This of course implies that δG µν gravitates as a massless graviton and is consistent with this being the massless eigenstate since the Einstein-Hilbert structure of GR is fixed uniquely for such a field. It is also consistent with our general arguments of section 4.1.1.
The second column shows a possible direct decay channel for the massive spin-2 field into two gravitons. Remarkably, these interactions are absent and therefore such a decay is not possible. Again this is consistent with the general arguments of the previous section. Of course, there may still be graviton production due to decay of the massive spin-2 field mediated by SM interactions, but these will generically be heavily suppressed.
The third column displays the gravitational interaction between the massive and massless spin-2 fields. It also captures the tree-level process of the inverse decay discussed later on with respect to possible production mechanisms for the massive spin-2 field. We note that these terms have no α dependence, which already indicates that the massive spin-2 field gravitates with the same strength as SM particles. That the gravitational stress-energy tensor of the spin-2 field indeed coincides with the one obtained via the Noether procedure in flat space follows directly from the general results of ref. [62] which we review in detail in appendix C. For a confirmation of these arguments through an explicit calculation, see [19].
Finally, the last column displays the self-interactions of the massive spin-2 field. Here we note that there are a variety of terms present and all come with factors of α. In particular, in the small α limit, some of these self-interactions will be enhanced as compared to standard GR. This enhancement is particularly strong when α is small and m FP is large. Table 2. Coefficients of quartic interaction vertices (numerical factors neglected) in units of m −2 Pl . Vertices with a dimensionless coefficient are associated with two derivatives. In this table we have used the notation β = α 2 m 2 Pl (β 1 + β 2 )/(1 + α 2 ) and m 2 g = m 2 Pl /(1 + α 2 ).

DM phenomenology
We have demonstrated that the heavy spin-2 field of bimetric theory acts as a perfect DM candidate: it interacts extremely weakly with SM particles and gravitates in the same way as ordinary matter does. In addition to these features, for the model to be viable, we need to make sure that the heavy spin-2 field is stable, at least on cosmological time scales, and that its relic abundance can match the one inferred from current observations.

Production mechanisms
We first discuss production mechanisms for the heavy spin-2 field abundance possibly active in the early Universe. Given the extremely weak interactions between the heavy spin-2 field and SM matter, the usual scenario in which the DM relic abundance is built via the freeze-out mechanism cannot be realised since the heavy spin-2 particle is never in thermal equilibrium in the early Universe. This can be straightforwardly seen by comparing the Hubble rate H to the interaction rate Γ ∼ n DM /m 2 Pl , where n DM ∼ ρ DM /m FP ∼ Ω DM H 2 m 2 Pl /m FP is the DM number density (energy density, density fraction), with the Hubble parameter itself. The DM abundance is fixed by observations to Ω DM 1 and this means that for m FP ≥ H thermal equilibrium could never be realised. In other words, the Hubble stretch always dominated over the relevant interaction rate H Γ and the current DM abundance could not arise via the freeze-out of a thermal DM population.
Gravitational Production. One possibility is that the heavy spin-2 field could be efficiently produced at the end of inflation, due to the non-adiabaticity of the transition between inflation and the hot Universe. This mechanism is known as gravitational particle production, see [63][64][65], and is quite independent of the details of this transition. We will briefly explain why this scenario is not consistent with our setup. Gravitational production is most efficient for masses on the order of (or actually, slightly larger than) the Hubble parameter at the end of inflation H e . The DM abundance today can be written schematically as, where h ∼ 0.7 is the little Hubble parameter and T rh the reheating temperature. Notice however that a large out-of equilibrium component of heavy DM will generate isocurvature perturbations, which are strongly constrained by CMB data [66]. The bound on isocurvature perturbations translates into a lower limit on m FP , or alternatively, on the scale of inflation at thermalisation H e : m FP /H e 5. This implies that, in order to generate a DM abundance satisfying Ω DM h 2 ∼ 0.1, we find, m FP ∼ 10 14 10 7 GeV T rh 1/2 GeV . Freeze-in mediated via s-channel exchange of the massive spin-2 δM ; the α factors cancel out and the amplitude is the same as for the massless δG mediation. being the ratio of tensor-to-scalar primordial perturbation amplitudes measured by the Planck experiment [66]. Realistically, the reheating temperature will be at least a factor of a few lower, and the heavy spin-2 mass can be estimated to be at least as heavy as 10 10 GeV. As we will show in the following, this requires α 1, otherwise the spin-2 particle would decay too rapidly. Recall, however, that we cannot take α to be arbitrarily small because our perturbative expansion would otherwise break down as discussed in detail in Sec. 4.1.2. In fact, we find that this precludes the possibility of gravitational DM production within our framework: As shown in Fig. 4, the region in the (α, m FP ) parameter space where gravitational production is successful is actually excluded by our perturbativity condition.
Freeze-in. Even if thermal equilibrium is never attained, it is possible to populate the Universe with a nearly decoupled species via a slow "leakage" of the thermal bath. This is the so-called freeze-in mechanism [67], which results in a non-thermalised sector, composed of the heavy spin-2 particles in our case. In this setup two SM particles from the thermal bath annihilate and produce a heavy spin-2 pair via s-channel graviton exchange. Depending on the dynamics of reheating, the generation of DM can proceed either during reheating or in the following radiation dominated era, see [68,69]. This process is very slow and never counterbalanced by the opposite reaction because the heavy spin-2 abundance remains well below the thermal one at all times.
In our setup, in addition to the usual massless graviton δG exchange channel, freeze-in can also proceed via exchange of the heavy spin-2 field δM itself. The two channels give identical results since the α suppression for the SM SM → δM vertex is compensated by the 1/α enhancement of the δM self-interaction δM 3 . These production channels are illustrated in Figure 2.
The generation of DM can be described by a system of coupled Boltzmann equations as in [70], where the thermally averaged cross section is given by σv ∼ T 2 /m 4 Pl [68,69]. The only relevant difference between the two possible production epochs (reheating or radiation domination) is the scaling of the Hubble rate H ∼ ρ 1/2 : in the first case ρ ∝ a −3/2 whereas in the second case ρ ∝ T 2 ∝ a −2 . Depending on the efficiency of the reheating process, generally parametrised as 2 rh = π 2 g * T 4 rh /90m 2 Pl H 2 e ≤ 1 with g * = 106.75 being the total number of relativistic degrees of freedom during reheating (which we take to be those of the SM only), the ranges of spin-2 masses for which it is possible to generate the correct DM abundance are, One can also estimate the total DM abundance directly in radiation domination: matching the observed DM abundance Ω DM via freeze-in in this case means [69], where m p is the proton mass, Ω b the abundance of baryons, and η b ≈ 10 −9 the baryon asymmetry. Once again, since the scale of inflation cannot be too high in order to avoid overproduction of tensor modes (not observed in the CMB), this implies that the heavy spin-2 mass will be constrained to the range, 1 TeV m FP 10 11 GeV . (5.5) This shows that, in principle, freeze-in is a possible production mechanism for our spin-2 DM. However, we still need to combine this α-independent result with the requirement of perturbativity. As an illustration of the discussion at the end of Sec. 4.1.2, consider the following vertices and their contribution to the production via freeze-in, cubic: It is clear that, even for energies below the perturbativity bound, E < α m Pl , both the δM 3 vertex and the δM 4 vertex dominate over the cubic one, δGδM 2 , for α < 1. In order to ensure that the above cubic terms are dominant over the quartic ones, we would have to impose the stronger bound E < α 2 m Pl . Were we to go further in the expansion, demanding a "safe" order α suppression for everything beyond the cubic order we would recover the condition E < α 7 m Pl as derived in section 4.1.2. However, when computing the actual production amplitudes, we need to look at the complete diagrams, SM SM → δG → δM δM and SM SM → δM → δM δM , and compare them to SM SM → δM → δM δM δM . Then, as we already mentioned, the contribution of the δM 3 vertex to the production rate is the same as the δGδM 2 one. Moreover, the diagram with the quartic vertex has an additional factor of E α m Pl . Hence our perturbative bound E < α m Pl is enough to trust our expressions derived with the cubic vertices only. Higher order vertices will only contribute a finite number of corrections to our estimates, proportional to increasing powers of E α m Pl < 1, and can thus be safely ignored.

Decay and possible signatures
Since the heavy spin-2 particle does not carry any of the SM charges (which automatically follows from the blindness of gravity to said quantum numbers), it decays universally into all the kinematically allowed channels, i.e. into all SM particles X with masses m X ≤ m FP /2. The universality of the decay processes is a feature that our bimetric DM model shares with, for instance, Kaluza-Klein DM [71]. However, interestingly, in bimetric DM the massive eigenmode δM µν cannot decay into massless modes. In other words, there is no graviton production, nor gravitational waves signals, associated to bimetric DM decay.
The most obvious upper bound on the mass m FP comes from imposing that the DM be stable on cosmological timescales. Requiring that its lifetime exceeds the age of the Universe τ U = 13.8 Gyr implies α 2/3 m FP 0.13 GeV. From this constraint we can then derive a consistency upper bound on the DM mass within our perturbative framework. Our expansion (as well as the expression for the width in eq. (5.6)) is valid for m FP ≤ α m Pl . This limit intersects the bound on the lifetime at m FP ≈ 6.6 × 10 6 GeV. Consequently, as remarked before, gravitational particle production is not a viable mechanism to generate the DM as it only operates efficiently for much higher masses. Furthermore, the viable range for production via freeze-in is shrunk to, 1 TeV m FP 6.6 × 10 3 TeV . (5.11) Given this mass range of the heavy spin-2 field, we can search for distinguishing indirect decay signals. In fact, even tighter constraints than that of eq. (5.11) can be derived by using the (non)observation of SM particle fluxes due to DM decay in different channels. In general, the constraints on the individual decay widths are heavily dependent on the mass and the propagation properties of the primary and secondary decay products, see for example [72]. We gathered the most stringent constraints from DM indirect detection experiments in Fig. 3, where we show the bounds on the inverse partial decay widths as a function of the DM mass. At low DM mass, the strongest constraints for our model come from the Fermi LAT searches for γ-ray lines [5]: these are the strongest constraints overall, hovering over the 10 29 s for the DM lifetime, but only apply up to masses of the order of a TeV. In the intermediate region, for which TeV m FP 10 TeV, the most competitive limits come from the antiproton flux measured by PAMELA [73] instead; the fluxes obtained by the AMS-02 experiment are in the same range [74]. Moreover, the constraints from the Extragalactic Gamma Ray Background of DM decaying into all SM channels are also in the same ballpark, see [75] we report here only the most significant ones from the muonic, tauonic, and bottom quark channels. Finally, for the highest mass range we are interested in, m 10 TeV, the searches for neutrino lines in IceCube provide the most relevant limits, around 1/Γ νe 10 28 s [76].
Roughly speaking, we can see that the limit obtained for a DM mass in the range (5.11) is approximately 10 orders of magnitude stronger than the bare limit coming from the lifetime of the Universe. In the perturbative regime this translates into an upper limit on the mass of 1 TeV m FP 66 TeV . (5.12) This very limited mass range for heavy spin-2 DM is one of the predictions of our model: a measured DM mass within this narrow range would be a strong indication in support of this model. In Fig. 4 we collect the strongest constraints on the total decay width, mostly coming from the DM decay into photons and neutrinos, together with the perturbativity limit. We check them against the different mass ranges available for freeze-in production. The plot shows the available (α, m FP ) parameter space for bimetric spin-2 DM: the mass of the heavy spin-2 particle is constrained to be in the 1 to 100 TeV range, while the value of α is approximately between 10 −11 and 10 −15 . Translating the latter to the value pertaining to the massive spin-2 self-interactions, we find this scale to be in the 10 3 GeV to 10 8 GeV range.
To summarise: The highest and lowest viable values for the mixing parameter α come from intersecting the available mass range for freeze-in production in (5.5)  x /s EGRB, µμ EGRB, ττ EGRB, bb PAMELA p, bb or tt AMS p/p, bb Fermi LAT γ-ray lines IceCube, ν e ν e Figure 3. Constraints on the partial decay width from analysis of the EGRB for DM decaying to muons (dashed yellow line), taus (solid red) and b quarks (dashed blue) [75]. The purple solid line and the violet dashed one show the bounds imposed by the antiproton measurements of the PAMELA [73] and AMS-02 experiments [74], respectively. We indicate with a solid green line the constraint due to Fermi LAT searches for γ-ray lines [5] and, with a black dashed line, the bound due to the observation of the high-energy electronic neutrino flux by the IceCube experiment [76].
limit on the DM mass m FP is also obtained by joining the latter two constraints: the DM must be stable enough and the energy at production must be within perturbativity.

Discussion
A heavy massive spin-2 field, whose gravitational interactions are described by the ghost-free bimetric theory, possesses all the desired features of a DM candidate. The extremely weak coupling of the spin-2 field to SM matter furthermore explains the absence of DM signals in dedicated detection experiments and collider searches. The Planck mass αm g ∼ αm Pl of the second metric can range from 1 to 10 4 TeV and the Fierz-Pauli mass for the spin-2 field is on the order of m FP ∈ [1,66] TeV. This narrow mass range for the DM candidate is one of the most distinct features of our model. Note however that the upper bound was obtained from the requirement of remaining in the perturbative framework. In principle, non-perturbative methods (which are presently unknown) could reveal that a larger spin-2 mass is also consistent with phenomenology.
Another exceptional property is the universal decay of DM into all SM particles along with the absence of a decay channel into massless gravitons at tree level. Our scenario predicts that mass and interaction scale of the heavy spin-2 field are of the same order of magnitude and only slightly larger than the weak scale. The largeness of the physical Planck mass m Pl is responsible both for suppressing the interactions of DM with baryonic matter and for bringing the theory close to GR. We have therefore not created any new hierarchies of energy scales and moreover related the puzzle of a large Planck scale to the extremely weak interactions of DM.
The above constraints on bimetric parameters are consistent with all current gravity tests. They correspond precisely to the overlap of the regions m 2 FP Λ for the spin-2 mass and α 1 for the ratio of Planck masses discussed in Sec. 3.2. Indeed, the Compton wavelength of the spin-2 field is tiny, approximately between 10 −19 cm and 10 −17 cm, resulting in typical Vainshtein radii of r V ∼ 10 −10 cm for the Sun, and r V ∼ 10 −24 cm for millimetre/submillimetre tests of the gravitational inverse-square law. These values ensure the validity of the linear approximation for all possible local gravity tests [77,78]. Relative corrections to GR solutions (the ratio of the fifth force to the Newtonian force) involve a factor of α 2 exp(−m FP r), which is extremely small in our setup. The largest relative deviation from Newton's law, which would be accessible through the sub-millimetre tests, is thus at most ∼ 10 −22 exp(−10 16 ), which is far beyond the reach of any experiment. In a very similar fashion, the smallness of the parameter α ensures that bimetric predictions for cosmology are essentially indistinguishable from the ΛCDM model and that perturbations remain stable. Bimetric theory in the above parameter region therefore resembles GR, but with an additional tensor field behaving like cold DM.
The fact that static spherically symmetric solutions in the weak field approximation resemble GR solutions for α 1 and m 2 FP Λ may suggest that in this parameter range black holes do not have any specific features which would distinguish them from those of GR. For instance, the instability found for spherically symmetric black holes [79,80] -which is clearly a distinctive feature of bimetric theory, since in GR such black holes are stabledisappears for large m FP , since the instability range is limited to m FP r −1 S . On the other hand, the absence of Birkhoff's theorem in bimetric theory allows for the existence of hairy black holes, see [81,82] for particular examples and [83] for a recent discussion. It is therefore feasible that hairy black holes are present in our scenario, possibly giving rise to distinct observational features of bimetric theory. It remains to be answered how different these hairy black holes can be from GR black holes, but this question lies beyond the scope of our paper.
Another distinct property of spin-2 DM are its enhanced self-interaction terms whose form is fixed by the ghost-free structure of the bimetric potential. Self-interacting DM is known to produce observable effects in collisions of galaxy clusters but so far the constraints are not very stringent [84]. Moreover, DM self-interactions would induce distortions in the DM power spectrum at small scales while leaving the baryonic one untouched. However, these interactions are mediated by δM itself, which confines the effectiveness of the corresponding force to risible length scales in astronomical or cosmological terms. Thus, current constraints on the self-interaction cross section are of little relevance in our case, despite the strength of the 1/α enhancement, since in practice the DM particles never feel each other.
Let us finally point out that, in addition to the aforementioned characteristics of our model, the decay of spin-2 DM would in principle exhibit a slightly different spectrum of secondaries (especially neutrinos), compared to the decay of, for instance, a scalar singlet. Indeed, since the δM µν field couples directly to the energy-momentum tensor of the SM, the two-vector final state (such as ZZ or W + W − ) will be mostly transversal, that is, will carry spin 2. This is in contrast to DM being a scalar, since in that case only longitudinal states can be produced (at tree level, but other polarisations can appear through higher dimensional operators). It can also differ from models with a Kaluza-Klein massive graviton, since in those constructions the transversal-to-longitudinal ratio can be different (depending on the structure of the extra dimensional setup), see for example [85]. However, this is a very model-dependent statement and we emphasise that in bimetric theory all predictions are fixed by demanding consistency of the theory. Once the vector bosons decay, the final spectra of secondary neutrinos will bear information about the spin in the guise of peculiar spectral features, see for instance the discussion in [86]. These features, however, are very small and hardly within the reach of current experiments. This theory suffers from the vDVZ discontinuity: in the limit m FP → 0 the solutions to the equations following from (A.2) do not correspond to those of linearised GR. The discontinuity manifests itself in a solution with a point-like source. For an ansatz for the metric perturbation of the form (3.1), the equations of motion derived from (A.2) imply, where C 3 is an integration constant, to be fixed by the matching to the source. The same result can be derived from (3.3) by applying the limit (A.1). For r m −1 FP the gravitational force is exponentially suppressed by the Yukawa potential. At smaller radii, r m −1 FP , the metric functions exhibit the correct r −1 power-law behaviour, but the ratio of the metric functions is |ν/λ| = 2, whereas in GR this ratio is exactly unity. This deviation persists for all radii satisfying r m −1 FP and thus constitutes a discontinuity of the zero-mass limit. As a consequence, the Fierz-Pauli theory does not pass the most basic Solar System tests which constrain |ν/λ| to be close to unity.

B The cubic action
Here we provide the explicit form of the bimetric action expanded up to cubic order in fluctuations around the proportional backgrounds.

B.1 Metric fluctuations
We consider fluctuations around the vacuum solutionf µν =ḡ µν , in the form, The metric determinant expanded up to cubic order is given by, where the square brackets around a tensor denote the trace, e.g.
[h] ≡ḡ µν h µν etc. Obviously the analogous expression for |f | is obtained by the formal replacement h → . We also write down the expansion of the square root matrix S = g −1 f , which on the proportional backgrounds can be obtained from a formal expansion of the form, where the binomial coefficient is given by 1/2 n = (1/2) n /n! (with the Pochhammer symbol (1/2) n representing a falling factorial). Written out in full the expansion to cubic order is given by, Recall the bimetric action absent matter sources, After a lengthy exercise in algebra we find that the complete cubic action, modulo boundary terms, can be written as, Here the first term is just a constant piece which is present for constant curvature backgrounds but does not contribute to the equations of motion. The remaining terms in the first line are given by the following GR-like terms, originating from |g|(R(g) − 2Λ) + α 2 |f |(R(f ) − 2Λ) with Λ = α 2 m 2 g (β 0 + 3β 1 + 3β 2 + β 3 ) = m 2 g (β 4 + 3β 3 + 3β 2 + β 1 ), L The second line in (B.6) contains terms from the interactions. First the mass term is, L where β ≡ α 2 m 2 g (β 1 + β 2 ).

C Stress energy tensors
Here we provide a general argument for the flat space on-shell physical equivalence between Noether and gravitational stress energy tensors, following [62]. Consider a gravitational Lagrangian of the following general form, L Tot = L GR (G) + L m (G, ∂G, δM, ∇δM ) . (C.1) As our notation suggests, the total Lagrangian L Tot consist of a gravitational part L GR for the nonlinear metric G which has the form of the standard Einstein-Hilbert term and a "matter" part L m . The matter field δM couples to the metric G as well as its first derivatives (which must enter in the form of Christoffel symbols due to covariance). As we discussed in section 4.1.1, the above Lagrangian captures the structure of bimetric theory expanded to infinite order around proportional backgrounds in terms of the linear massive mode δM and the nonlinear gravitational metric G µν =ḡ µν + δG µν /m Pl . The general results derived in the following therefore all hold true in our setup.
We compute the gravitational stress energy tensor corresponding to the above theory from, Note that the last term is usually absent since standard matter only couples to the graviton and not its derivatives.
On the other hand, the Noether stress-energy tensor derived from translational symmetry in flat space is, As a final ingredient we recall the Euler-Lagrange equations for the matter field in the form, Now we evaluate the gradient, Using the definition (C.2) and the equations of motion (C.4) we can use this to obtain, Further manipulation and using the covariant conversation law ∇ ρ T ρ ν = 0 finally gives, This is the general correspondence between the stress energy tensors in curved spacetime for a theory of the form we are considering. Note that in general the stress energy tensors are of course not equal and furthermore, the Noether stress energy tensor is strictly speaking derived from flat space considerations. In order to check the physical implications of (C.7) it is instructive to study the momenta defined by, We see that (C.7) implies that (the constant of integration can be shown to vanish), Therefore, in flat space where ∂ ν G µσ = 0, we find that π ν = P ν which shows the physical equivalence between the prescriptions. From (C.7) it also follows that in flat space (and/or for theories where matter does not couple to the derivative of the graviton) we always have (again the constant of integration can be shown to vanish), Hence, in flat space the integrated stress energy tensors are always equal on-shell (modulo a factor of 2 which is due to our definition of T µν ). This last relation has an immediate and interesting implication. Namely, considering a non-relativistic gas (fluid) of matter satisfying Bose statistics, the Noether stress energy is known to correspond to a stress energy of the form of dust τ ρ ν ∼ diag(ρ, 0, 0, 0). The above relations now imply that this is also the form of the source for the gravitational field and therefore the massive spin-2 field will behave just as cold DM non-relativistically. This can also be explicitly verified directly from our expressions for the cubic interaction terms in appendix B. 10