# Recent developments in bimetric theory

Angnis Schmidt-May, Mikael von Strauss

► To cite this version:

Angnis Schmidt-May, Mikael von Strauss. Recent developments in bimetric theory. Journal of Physics A: Mathematical and Theoretical, 2016, 49, 10.1088/1751-8113/49/18/183001. [https://doi.org/10.1088/1751-8113/49/18/183001](https://doi.org/10.1088/1751-8113/49/18/183001). [arXiv:1605.08457](https://arxiv.org/abs/1605.08457)

HAL Id: hal-03645167

https://hal.science/hal-03645167

Submitted on 23 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
TOPICAL REVIEW • OPEN ACCESS

Recent developments in bimetric theory

To cite this article: Angnis Schmidt-May and Mikael von Strauss 2016 *J. Phys. A: Math. Theor.* **49** 183001

View the [article online](#) for updates and enhancements.

You may also like

- [The mean gauges in bimetric relativity](#)
  Francesco Torsello

- [Structures in multiple spin-2 interactions](#)
  Oliver Baldacchino and Angnis Schmidt-May

- [Initial data and first evolutions of dust clouds in bimetric relativity](#)
  Mikica Kocic, Francesco Torsello, Marcus Högás et al.

<figure><img src="image_2.png" /></figure>

IOP ebooks™

Bringing together innovative digital publishing with
leading authors from the global scientific community.

Start exploring the collection-download the
first chapter of every title for free.
## Topical Review

# Recent developments in bimetric theory

Angnis Schmidt-May<fnref n="1" />,<fnref n="3" /> and Mikael von Strauss<fnref n="2" />

E-mail: angniss@itp.phys.ethz.ch and strauss@iap.fr

Received 10 December 2015, revised 2 March 2016

Accepted for publication 7 March 2016

Published 29 March 2016

<figure><img src="image_3.png" /></figure>

CrossMark

### Abstract

This review is dedicated to recent progress in the field of classical, interacting, massive spin-2 theories, with a focus on ghost-free bimetric theory. We will outline its history and its development as a nontrivial extension and generalisation of nonlinear massive gravity. We present a detailed discussion of the consistency proofs of both theories, before we review Einstein solutions to the bimetric equations of motion in vacuum as well as the resulting mass spectrum. We introduce couplings to matter and then discuss the general relativity and massive gravity limits of bimetric theory, which correspond to decoupling the massive or the massless spin-2 field from the matter sector, respectively. More general classical solutions are reviewed and the present status of bimetric cosmology is summarised. An interesting corner in the bimetric parameter space which could potentially give rise to a nonlinear theory for partially massless spin-2 fields is also discussed. Relations to higher-curvature theories of gravity are explained and finally we give an overview of possible extensions of the theory and review its formulation in terms of vielbeins.

Keywords: interacting spin-2 fields, massive gravity, bimetric theory, modified gravity

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 3.0 licence](https://creativecommons.org/licenses/by/3.0/). Any further distribution of this work must

maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Institut für Theoretische Physik, Eidgenössische Technische Hochschule Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland</footnote>

<footnote n="2">² UPMC-CNRS, UMR7095, Institut d'Astrophysique de Paris, GReCO, 98bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3">³ Author to whom any correspondence should be addressed.</footnote>
# 1. Introduction

## 1.1. Motivation

The Standard Model of particle physics contains massive and massless fields with spin 0, 1/2 and 1. Gravitational interactions are attributed to a spin-2 field which, in the standard framework of general relativity (GR), is massless and possesses nonlinear self-interactions. Even though the Standard Model and GR are experimentally and observationally well-tested, several phenomena still remain unexplained and motivate the study of theories beyond the standard picture. In particular, two of the biggest unresolved problems concern the presently unknown nature of dark matter and dark energy and, in order to account for these constituents in a satisfactory way, introducing new physics becomes unavoidable. Additional degrees of freedom could be of the same type as the fields already present in the standard scenarios or, more interestingly, they could arise from heretofore unknown field theories. While the field theories for spin 0, 1/2 and 1 are well-understood, the treatment of higher spins turns out to be much more difficult. One of the simplest, or at least most natural, new ingredients that could be added to the known models is a massive spin-2 field whose presence is expected to mostly affect the gravitational sector. This may be desirable since modifying gravity is motivated by the fact that the Standard Model of particle physics is based on the very solid framework of a renormalisable quantum field theory, while a quantum theory of gravity does not yet exist and hence GR is not expected to be complete. Moreover, both the dark energy and the dark matter problems are intimately related to gravity but cannot be solved in the context of GR without raising additional questions.

Interactions for massive spin-2 fields have long been thought to inevitably give rise to ghost instabilities and only recently a ghost-free theory for nonlinear interactions between massive and massless spin-2 fields has been found. Since the construction of nonlinear massive gravity and its extension to bimetric theory there has been significant progress in the field, both on the theoretical and the phenomenological side. Two long review articles on the subject have already been written by Hinterbichler [1] and de Rham [2]. While these references focus on massive gravity, its gauge invariant formulation in terms of Stückelberg fields and its cosmology, this review is mostly dedicated to the manifestly covariant and dynamical bimetric theory. We focus on theoretical aspects and the structure of the theory (which is enforced by consistency) and mention its application to cosmology only as an aside. Moreover, all our considerations in this article will be at the classical level, although quantum corrections have been under preliminary investigations and are discussed in more detail in de Rham's review [2], mostly in the context of massive gravity.

The following subsection summarises the interesting historical developments that led to the construction of consistent spin-2 interactions. It is no prerequisite for the subsequent discussions and the reader more interested in the details of massive spin-2 theories may directly proceed to section 2.

## 1.2. Historical background

Since massive gravity and bimetric theory have a rather long history starting with the work by Fierz and Pauli in 1939, it is impossible to give sufficient credit to all the groups that have contributed to the field over the last 75 years. We try here to collect the most important inputs towards the modern construction and pay attention to historical accuracy to the best of our knowledge.
1.2.1. Early attempts. The program of investigating massive spin-2 fields was initiated by Fierz and Pauli, who derived the unique classically consistent linearised theory of a free massive spin-2 field $h_{\mu\nu}$ propagating in Minkowski space–time [3, 4]. They demonstrated that the corresponding Lagrangian is of the form

$$ \mathcal{L}_{\text{FP}} = \frac{1}{4}(-\partial_{\mu}h_{\alpha\beta}\partial^{\mu}h^{\alpha\beta} - 2\partial_{\mu}h\partial_{\nu}h^{\mu\nu} + 2\partial_{\mu}h^{\mu\nu}\partial^{\alpha}h_{\alpha\nu} + \partial_{\mu}h\partial^{\mu}h) - \frac{m_{\text{FP}}^2}{4}(h_{\mu\nu}h^{\mu\nu} - h^2), \quad (1.1) $$

where $h = \eta^{\mu\nu}h_{\mu\nu}$ and $m_{\text{FP}}$ is the mass parameter. The first line here corresponds exactly to the kinetic operator obtained by linearising the Einstein–Hilbert action of GR around flat space, while the second line encodes the quadratic non-derivative self-interactions of $h_{\mu\nu}$ that render it massive. The equations of motion derived from the above Lagrangian are equivalent to the system of equations

$$ (\Box - m_{\text{FP}}^2)h_{\mu\nu} = 0, \quad \partial^{\mu}h_{\mu\nu} = 0, \quad h = 0. \quad (1.2) $$

The first of these is a massive wave equation, while the latter two are constraints on components of $h_{\mu\nu}$. In particular, Fierz and Pauli showed that if the relative coefficient $(-1)$ of the two parts in the mass term of (1.1) is changed in any way, the on-shell condition of tracelessness is lost and a ghost-like scalar mode inside $h_{\mu\nu}$ becomes propagating. At the classical level a ghost is a field with negative kinetic energy which gives rise to an unbounded Hamiltonian and thus causes fatal instabilities; at the quantum level ghosts must be avoided in order to ensure unitarity. It is therefore crucial to work with the above Lagrangian with correct relative coefficient in the mass term. In four space–time dimensions, it describes the on-shell propagation of a traceless, transverse and symmetric tensor field $h_{\mu\nu}$ with five massive degrees of freedom. This allows us to identify $h_{\mu\nu}$ with a massive spin-2 field with helicities $\pm 2, \pm 1, 0$. In all that follows, when using the terms mass and spin, we refer to their correspondence to degrees of freedom of relativistic field equations; we will rarely speak about the quantum nature of these concepts. In a similar fashion, we will frequently employ the particle perspective when dealing conceptually with nonlinear generalisations of the Fierz–Pauli Lagrangian, such as bimetric theory.

During the following years much of the efforts in field theory was directed elsewhere, developing a firm understanding of GR and exploring the booming realm of particle physics. Not much consideration was paid to Fierz and Pauli’s theory of massive spin-2 until the early 1970s. We note however that related questions concerning bimetric theories were raised early on by Rosen [5, 6] and later by Aragone and Deser in [7] and by Isham, Salaam and Strathdee in [8]. A major development came about when van Dam and Veltman [9] and Zakharov [10] (see also [11]) independently investigated particular consequences of the Fierz–Pauli Lagrangian interpreted as a theory of a massive graviton. They realised that, in the presence of matter sources, the zero-mass limit of the theory is discontinuous, a property which is now referred to as the van Dam, Veltman and Zakharov (vDVZ) discontinuity. More precisely, this limit does not result in a theory for a single massless spin-2 field like linearised GR, but also contains a propagating scalar field which couples to the trace of the stress-energy source. The observational consequence would be an inferred difference in the bending of light around massive sources which was so severe that the theory would have been ruled out already at this time. Shortly thereafter, however, Vainshtein recognised a loophole in this reasoning [12]. He argued that, due to the presence of more scales in the theory when coupled to a source, a scalar mode becomes strongly coupled below some distance $r_V$ (the Vainshtein radius). Hence the linearised analysis breaks down in this regime and a nonlinear completion of the theory is necessary in order to address any questions at distances within $r_V$ in a consistent
manner. In particular these findings made it possible again to recover GR at short distances in the zero-mass limit. That this recovery can indeed be realised at the nonlinear level was demonstrated in [13] for the case of ghost-free bimetric theory (for a recent review on the Vainshtein mechanism see [14]).

As an immediate response to Vainshtein's idea, Boulware and Deser studied the consistency of a wide class of possible nonlinear extensions of the massive Fierz-Pauli theory. They concluded that it was inevitable to introduce an extra propagating ghost-like scalar mode in any nonlinear extension of the theory [15]. In fact, this scalar mode, the Boulware-Deser ghost, is the same mode that at the linear level was removed by the trace constraint on $h_{\mu\nu}$. The implication of this analysis was that no consistent nonlinear theory of massive spin-2 fields could exist. As we now know this strong conclusion was incorrect for two reasons: (1) Boulware and Deser did not consider the most general nonlinear extensions possible in their analysis. They considered only non-derivative self-interactions of $h_{\mu\nu}$ that were given through a general analytic function $f(h^2 - h_{\mu\nu}^2)$, which naturally Taylor-expands to the Fierz-Pauli mass term at lowest order. As we will see later, in the consistent theory the self-interactions are contained in very specific scalar functions (the elementary symmetric polynomials) constructed out of the matrix argument [$\sqrt{1 + \eta^{-1}h}]^{\mu}_{\nu}$. The correct field dependence of the interactions is therefore not of the form assumed in the proof by Boulware and Deser. (2) In their Hamiltonian analysis, they expected one specific equation (the equation of motion for the lapse variable of the metric) to provide the constraint that removes the ghost. The analysis of the consistent theory reveals that, in fact, it is a different equation (a rather contrived combination of lapse and shift equations) that gives the constraint.

The conclusions of Boulware and Deser's analysis were so widely accepted that no further progress was made in the field for another 30 years and it would be almost 40 years until a consistent theory was fully developed. In hindsight this is unfortunate since, in the meantime, some interesting ideas did not receive the deserved attention due to the strong no-go theorem. For example, the correct structure of interactions for massive spin-2 fields was in fact partly suggested very early on by Zumino and also Chamseddine [16, 17], but without addressing the ghost problem (in fact it can be noted that [16] actually predated the no-go theorem). Similarly, the correct structure in the vielbein formulation was partly written down in [18], again failing to address the ghost issue. It has also recently been pointed out in [19] that some attempts to construct a theory of nonlinear massive spin-2 fields were made by Maheshwari [20] using ideas from the works of Ogievetsky and Polubarinov [21, 22]. These constructions however went largely unnoticed and did not contribute towards the modern ideas and understanding of nonlinear massive gravity.

Another approach that has become popular in recent years and deserves mentioning is the attempt to construct consistent spin-2 mass terms that break Lorentz invariance [23–25]. Based on the experience gained from the recent progress made in constructing the Lorentz invariant theories, similar progress has also been made in constructing Lorentz-breaking theories of massive gravity [26–28]. In this review we will however restrict ourselves to the class of theories that respects Lorentz invariance.

1.2.2. A renewed interest. After the very precise confirmation of the accelerated expansion of the Universe in 1998 [29, 30], a huge effort was devoted towards understanding better the underlying physics of this discovery. Within the context of standard GR this observed accelerated expansion of spacetime itself requires the addition of an additional source term with rather strange behaviour as compared to standard matter sources. The main-stream philosophy at this time was to hold firm in GR and not worry too much about the origin of
such a source but simply collectively call this mysterious source ‘dark energy’. The simplest explanation for the acceleration is the presence of a constant source term in Einstein’s equations. Adopting this view led to the celebrated ‘Λ cold dark matter’ (ΛCDM) concordance model of cosmology, which is in excellent agreement with observational data thus far [31]. On the other hand, a growing community of theoretical physicists with a particle physics oriented mind was now realising the pressing nature of the cosmological constant problem [32]: the small value of the observed acceleration does not fit in with expectations from a particle physics perspective, where a cosmological constant is naturally associated with vacuum energy. Unless additional symmetries (such as supersymmetry) are at work, a natural value for the vacuum energy scale is the mass of the heaviest field in the theory, which in any scenario is many orders of magnitude higher than the observed value for the cosmological constant.

In addition to the poorly understood nature of dark energy, cosmologists seek to explain the presence of an unidentified matter component, commonly referred to as ‘dark matter’, which, in the context of GR, is required to account for the observational data at distances ranging from galactic to cosmological scales [31, 33, 34].

Since quantum field theory is such a rigorous framework and both the nature of dark matter and in particular dark energy seems so deeply connected with the large scale behaviour of gravity, which was really only tested within the solar system, many theorists started to look for an answer by modifying the gravitational sector of field theory.

In the beginning much attention was given to extra dimensional setups which were mainly inspired by the additional dimensions arising in string theory and by the braneworld scenarios geared towards addressing the Higgs hierarchy problem of the Standard Model as well as supersymmetry breaking through anomaly mediation [35–37]. With respect to the cosmological constant problem, a particular interest was paid to models of brane induced gravity [38–40] and similar constructions (e.g. cascading gravity [41, 42]). Brane induced gravity models, and in particular the codimension-one DGP model [38, 43] (for related investigations see [44]), were historically very important for a renewed interest in studying massive gravitons, since a generic feature of these models was the appearance of massive spin-2 resonance states on the brane that makes up our Universe in these models (see e.g. [45]). Due to the generic Yukawa suppression of the gravitational potential mediated by a massive field, a massive graviton also became interesting in itself for addressing the cosmological constant problem. The hope was that this exponential suppression could sufficiently weaken gravity at large distance scales to screen out a large vacuum energy coming from the matter sector and thus lead to a small effective cosmological constant. Even though this picture is correct in linear Fierz–Pauli theory, it eventually turned out that it would not work in the nonlinear theory without fine-tuning (as discussed in e.g. [46]).

In 2002 Arkani-Hamed, Georgi and Schwartz proposed a new perspective on studying effective theories with explicitly broken general covariance, in close analogy to symmetry breaking in spin-1 theories and the associated emergence of Goldstone modes [47]. These ideas were quite general but in particular provided a new language for analysing the internal consistency of massive gravity. More explicitly, the setup was based on intuition from the Goldstone boson equivalence theorem [48, 49], which relates the physics of longitudinal modes of spin-1 gauge bosons to the physics of Goldstone modes at high energies. The authors of [47] suggested a similar correspondence in gravitational theories and that, in certain energy regimes, the complex problem of self-interacting spin-2 fields could be simplified to studying only their scalar longitudinal components.

In 2005 Creminelli, Nicolis, Papucci and Trincherini followed up on the ideas of Arkani-Hamed, Georgi and Schwartz (further inspired by the results of [23, 24, 50]) and attempted to
explicitly construct a consistent theory of massive gravity using a bottom-up approach<fnref n="4" /> [51]. As a bitter irony of history (and humbling lesson in importance of rigor), despite a very beautiful analysis they reached the erroneous conclusion that the ghost problem of nonlinear massive gravity could not be resolved and found the Boulware–Deser ghost reappearing again. This result was based on an unfortunate sign mistake which arose from copying a basic equation of [47] (for which the sign was not important). It is worth mentioning that without this sign mistake the discovery of a consistent theory of massive gravity could have been made already in 2005.

### 1.2.3. Massive gravity rediscovered: de-Rham–Gabadadze–Tolley (dRGT) theory.
In 2009 the issue of massive gravitons was further pursued by Gabadadze who modified GR by introducing an auxiliary extra dimension [53]. This work was shortly followed up by (and with) de Rham in [54, 55]. Their approach was based on an interesting field theoretical tool to generate a mass term for a vector (or scalar) field by imposing boundary conditions in an auxiliary extra dimension. Even though introducing an extra unphysical dimension and fixing boundary conditions by hand seems rather ad hoc, the procedure itself straightforwardly extends to the spin-2 case. The massive spin-2 model obtained in this way was demonstrated to be ghost-free to cubic order in a ‘decoupling limit’ analysis [55]. Although the same setup was subsequently shown to be inconsistent in a fully nonlinear analysis by Hassan and Rosen [56], it provided important inspiration that pushed developments further. The belief that the auxiliary-dimension model was consistent motivated de Rham and Gabadadze to revisit the ghost analysis of nonlinear massive gravity by Creminelli et al.

In 2010 de Rham and Gabadadze studied generic extensions of the Fierz–Pauli Lagrangian (1.1) by higher-order interactions of the massive spin-2 fluctuation $h_{\mu\nu}$ [57]. Their analysis went to quintic order in the longitudinal component of the massive spin-2 field and demonstrated that its interactions could in fact be made ghost-free in a decoupling limit, correcting the conclusions of [51]. The decoupling limit analysis relies heavily on the aforementioned Goldstone boson analogy suggested by Arkani-Hamed, Georgi and Schwartz and requires taking a double scaling limit in order to study the dynamics of the longitudinal mode separately. As a follow up to [57], de Rham, Gabadadze and Tolley (henceforth dRGT) presented a nonlinear theory of massive gravity in whose decoupling limit they proved the absence of ghost for all nonlinear self-interactions of the longitudinal component [58]. The dRGT action is of the form

$$S_{\text{dRGT}} = m_g^2 \int d^4 x \sqrt{g} \left( R(g) + \frac{m^2}{2} \sum_{n=2}^{4} n! (4-n)! \alpha_n e_n(K) \right), \quad (1.3)$$

where the first term is the ordinary Einstein–Hilbert term of GR with Planck mass $m_g$ and the second term is the interaction potential for the graviton whose mass is set by the scale $m$. Furthermore, $\alpha_2 = 1$ while the two remaining $\alpha_n$ are arbitrary interaction parameters. The $e_n(K)$ are the elementary symmetric polynomials (see appendix A) constructed out of the<fnref n="5" /> matrix $K_\nu^\mu = \delta_\nu^\mu - [\sqrt{g^{-1}\eta}]_\nu^\mu$. Expanding the action in terms of $h_{\mu\nu} = m_g(g_{\mu\nu} - \eta_{\mu\nu})$ indeed results in a nonlinear extension of the Fierz–Pauli theory (1.1).

<footnote n="4">⁴ At the same time Deffayet and Rombouts independently used the same formalism to study the ghost and its relation to the Vainshtein mechanism [52].</footnote>

<footnote n="5">⁵ It should be noted that [58] did actually not present the action precisely in the above form. The interactions were resummed to give the square-root matrix structure, but it was not obvious from [58] that the sum in (1.3) would terminate at $n = 4$. This fact was pointed out by Hassan and Rosen in [46].</footnote>
The authors of [57, 58] also made the important observation of a loophole in the Hamiltonian analysis by Boulware and Deser [15] which had long been believed to forbid any consistent theory of nonlinear massive spin-2 interactions<fnref n="6" />. By avoiding the ghost in the interactions of the longitudinal component, the dRGT model satisfied one of the necessary requirements on the complete theory to be consistent and consequently became the most promising candidate for a consistent theory of massive gravity. However, as noted, the proof of absence of the Boulware–Deser ghost in [58] was not made for the full theory since the analysis focused on self-interactions of one scalar component only. One reason for the difficulty of conclusively proving the consistency was the complexity of the nonlinear spin-2 interactions whose structure required better understanding before a full Hamiltonian analysis could be performed.

### 1.2.4. Important generalisations.
The dRGT formulation of massive gravity was constructed with a perturbative expansion around a fixed flat background in mind. In 2011 Hassan and Rosen presented a reformulation of the dRGT action, which clarified the non-perturbative structure of the theory and identified consistent generalisations to more general backgrounds that were not apparent in the original formulation. Using basic properties of the elementary symmetric polynomials, they showed that the dRGT massive gravity action (4.7) can be reformulated and then generalised to [46]

$$S_{HR} = m_g^2 \int d^4 x \sqrt{g} \left( R(g) - 2m^2 \sum_{n=0}^{4} \beta_n e_n (\sqrt{g^{-1}f}) \right), \quad (1.4)$$

where $\beta_n$ are five arbitrary parameters whose role will be explained in detail later. Most notably, the new formulation involves an arbitrary (but fixed) reference metric $f_{\mu\nu}$, which in dRGT is strictly taken to be $\eta_{\mu\nu}$. It is worth to stress that even though the structures of (1.3) and (1.4) are very similar, there is no obvious way of getting to (1.4) from the original dRGT formulation in [58] (whereas obtaining (1.3) from (1.4) is straightforward) and thus the new formulation truly represents an important generalisation of the original massive gravity theory. In particular, the new formulation is indispensable for addressing any question about massive spin-2 interactions on a curved background. Even more importantly, as we will see, it suggests how to arrive at a fully dynamical theory of interacting spin-2 fields.

Going beyond the decoupling limit analysis of [58], Hassan and Rosen quickly utilised the formulation (1.4) to give a fully nonlinear consistency proof for dRGT massive gravity with flat reference metric [59]: in a Hamiltonian analysis based on the Arnowitt, Deser, and Misner (ADM) formalism [60], they showed that the complete nonlinear theory gives rise to a constraint that removes the Boulware–Deser ghost. Shortly after this, the proof was generalised to the case of an arbitrary reference metric $f_{\mu\nu}$ in [61]. These consistency proofs (see also [62, 63] and [27, 64–69] for subsequent confirmations) for the massive gravity action (1.4) essentially completed the program of finding a Lorentz-invariant theory for a massive spin-2 field initiated by Fierz and Pauli in 1939.

### 1.2.5. Hassan–Rosen bimetric theory.
A very important outcome of the generalised investigations of Hassan and Rosen laid in obtaining an extension of the massive gravity theory in which the reference metric $f_{\mu\nu}$ receives its own dynamics. As a consequence, the two metrics $g_{\mu\nu}$ and $f_{\mu\nu}$ are treated on the same footing in this bimetric theory and all fields in the action are determined dynamically. Bimetric theories of gravity have been subject of

<footnote n="6">⁶ In fact, to our knowledge, this caveat was first encountered in [53].</footnote>
earlier investigations in, for instance, [6–8, 17, 70–75] but, just as nonlinear massive gravity, they generically suffer from the Boulware–Deser ghost instability.

Shortly after the nonlinear theory for massive gravity had been developed, the unique ghost-free bimetric theory was presented by Hassan and Rosen in [76]. Its form is reminiscent of the massive gravity action

$$S_{HR} = m_g^2 \int d^4 x \sqrt{g} R(g) + m_f^2 \int d^4 x \sqrt{f} R(f) - 2m_f^4 \int d^4 x \sqrt{g} \sum_{n=0}^4 \beta_n e_n(\sqrt{g^{-1}f}), \quad (1.5)$$

but it includes an Einstein–Hilbert term of $f_{\mu\nu}$, in which $m_f$ is the ‘Planck mass’ for the second metric. Most importantly, bimetric theory provides a covariant formulation for massive spin-2 fields, in which now both metrics are dynamical and the structure of the action is in fact symmetric with respect to the interchange of $g_{\mu\nu}$ and $f_{\mu\nu}$. The consistency of bimetric theory was first demonstrated in [62, 76] and, as we shall see later, this established the first theory describing consistent nonlinear interactions of massive spin-2 fields with massless ones.

### 1.3. *Outline of the review*

We have chosen the following structure for this review: our starting point will be the linear theories for massless and massive spin-2 fields in flat background in section 2. The possibility to generalise the massless theory to arbitrary backgrounds by constructing the nonlinear theory of GR motivates us to look for a nonlinear completion of the mass term. In section 3, after reviewing the ghost problem, we pursue this goal by introducing a suitable set of ADM variables and provide a constructive proof for the consistent spin-2 interaction potential. Section 4 summarises the main features of ghost-free nonlinear massive gravity with general and flat reference metric, in both gauge fixed and in gauge invariant form. We then move on to the fully dynamical bimetric theory in section 5 where we write down its mass spectrum and also introduce couplings to matter. Classical solutions to the bimetric equations are discussed in section 6, first in general and then with a focus on black hole and cosmological solutions. In section 7 we review the phenomenon of partially massless (PM) spin-2 fields on de Sitter space and present the idea of realising partial masslessness at the nonlinear level. In this context, we also reveal a connection between bimetric theory and certain higher-derivative theories of gravity. Possible extensions of bimetric theory, including new kinetic terms in higher dimensions as well as multiple spin-2 interactions, are discussed in section 8. Finally, we conclude this review article with a list of several open questions in section 9.

## 2. Spin-2 fields on fixed backgrounds

Having outlined the history of massive gravity in the introduction, we now turn to the detailed description of spin-2 field theories. In this section we first discuss the linear theory for a massless spin-2 field in a flat background and see how the correct number of physical degrees of freedom emerges due to gauge invariance. The nonlinear completion of the corresponding Lagrangian is the Einstein–Hilbert action of GR from which one can derive the linear theory around general curved backgrounds. Thereafter we review the linear theory of massive spin-2 in flat space, where the presence of constraints lead to the correct number of propagating modes. The quest for a generalisation to curved backgrounds again requires a nonlinear completion of the spin-2 mass term and paves the way towards nonlinear massive gravity.
## 2.1. Massless spin-2 field

### 2.1.1. Flat space.
The Lagrangian for a massless spin-2 field $h_{\mu\nu}$ to quadratic order in the field in flat space with Minkowski metric $\eta_{\mu\nu}$ reads

$$ \mathcal{L}_{\text{lin}} = \frac{1}{2} h^{\mu\nu} \mathcal{E}_{\mu\nu}^{\rho\sigma} h_{\rho\sigma}, \qquad (2.1) $$

where the structure of the kinetic terms is captured by the two-derivative operator

$$ \mathcal{E}_{\mu\nu}^{\rho\sigma} = \frac{1}{2}(\eta_{\mu}^{\rho} \eta_{\nu}^{\sigma} \partial^{2} - \eta_{\nu}^{\sigma} \partial_{\mu} \partial^{\rho} - \eta_{\mu}^{\sigma} \partial_{\nu} \partial^{\rho} + \eta_{\mu\nu} \partial^{\sigma} \partial_{\mu} \partial_{\nu} - \eta_{\mu\nu} \eta^{\rho\sigma} \partial^{2}). \quad (2.2) $$

The corresponding equations of motion are

$$ \mathcal{E}_{\mu\nu}^{\rho\sigma} h_{\rho\sigma} = 0. \qquad (2.3) $$

As can easily be verified, the Lagrangian and equations for the massless spin-2 particle are invariant under the following gauge transformations

$$ h_{\mu\nu}(x) \mapsto h_{\mu\nu}(x) + \partial_{\mu}\xi_{\nu}(x) + \partial_{\nu}\xi_{\mu}(x), \qquad (2.4) $$

with vector gauge parameter $\xi_{\mu}(x)$. An equivalent manifestation of this gauge invariance is the linearised Bianchi identity: the divergence of the left-hand side of the equations of motion (2.3) is identically zero, i.e. $\partial^{\mu}\mathcal{E}_{\mu\nu}^{\rho\sigma}h_{\rho\sigma} = 0$.

As should be familiar from the spin-1 case (Maxwell's theory), we can use the symmetry transformation (2.4) to pick a convenient gauge. The de Donder gauge condition for the spin-2 field is $\partial^{\mu}h_{\mu\nu} - \frac{1}{2}\partial_{\nu}h = 0$. This is the analogue of the Lorenz gauge for massless vectors and constrains four of the ten components in $h_{\mu\nu}$. In this gauge the equations of motion assume the simple form

$$ \Box h_{\mu\nu} = 0. \qquad (2.5) $$

After fixing the de Donder gauge there still exist residual transformations, namely those with $\Box\xi_{\mu} = 0$, that leave the de Donder gauge intact. Since the four residual gauge parameters satisfy the same equation as the field, they can directly be invoked to remove further redundant degrees of freedom. In total, gauge invariance therefore removes eight of the ten components in the symmetric field $h_{\mu\nu}$, and the only propagating modes are the two degrees of freedom of a massless spin-2 particle.

It is also possible to couple the massless spin-2 field to other fields by introducing an external source $T_{\mu\nu}$ into the Lagrangian

$$ \mathcal{L}_{\text{lin}} = \frac{1}{2} (h^{\mu\nu} \mathcal{E}_{\mu\nu}^{\rho\sigma} h_{\rho\sigma} - \kappa h^{\mu\nu} T_{\mu\nu}), \qquad (2.6) $$

where $\kappa$ is a coupling constant of mass dimension minus one (here $h_{\mu\nu}$ is taken to have mass dimension one). The linearised Bianchi identity now implies $\partial^{\mu}T_{\mu\nu} = 0$, i.e. conservation of the source.

### 2.1.2. Curved space.
The graviton in GR is a massless spin-2 particle and should hence be described by the Lagrangian (2.6). On the other hand, the equivalence principle tells us that the field should couple to all kinds of energy in the same manner, including its own stress-energy tensor. Implementing this requirement in (2.6) iteratively introduces nonlinearities in $h_{\mu\nu}$ and eventually leads to the Einstein-Hilbert action for GR (see e.g. [77, 78])

$$ S_{\text{GR}} = M_{\text{Pl}}^2 \int d^4 x \sqrt{-g} R + \int d^4 x \sqrt{-g} \mathcal{L}_{\text{m}}, \qquad (2.7) $$
The corresponding equations of motion are Einstein's equations

$$R_{\mu\nu} - \frac{1}{2}g_{\mu\nu}R = \frac{1}{M_{\text{pl}}^2}T_{\mu\nu}, \qquad (2.8)$$

where $R_{\mu\nu}$ is the Ricci tensor with trace $R$ and the stress-energy tensor is derived from the matter Lagrangian $\mathcal{L}_m$ as

$$T_{\mu\nu} \equiv -\frac{1}{\sqrt{-g}} \frac{\partial(\sqrt{-g} \mathcal{L}_m)}{\partial g^{\mu\nu}}. \qquad (2.9)$$

The Einstein-Hilbert action is the gauge invariant nonlinear extension of (2.6). It can be defined uniquely, as the field theory that describes nonlinear self-interactions of a massless spin-2 particle. The gauge transformations (2.4) are linearised general coordinate transformations (GCTs), under which the metric at the nonlinear level transforms as

$$g_{\mu\nu} \mapsto g_{\mu\nu} + \nabla_{\mu}\xi_{\nu} + \nabla_{\nu}\xi_{\mu}. \qquad (2.10)$$

The presence of a gauge symmetry at the nonlinear level ensures that the full theory propagates the same number of degrees of freedom as its linear version.

The Einstein equations (2.8) allow for flat space solutions, $g_{\mu\nu} = \eta_{\mu\nu}$, in the case of vanishing matter source. Linearising the action around this background in the perturbation $h_{\mu\nu} = \kappa^{-1}(g_{\mu\nu} - \eta_{\mu\nu})$ results in the linear theory for a massless spin-2 field (2.6).

On the other hand, a non-vanishing background value of the source will give rise to a curved background metric $\bar{g}_{\mu\nu}$ and around this solution the action can be linearised as

$$\mathcal{L}_{\text{lin}}' = \frac{1}{2}(h^{\mu\nu}\bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma}h_{\rho\sigma} - \kappa h^{\mu\nu}\delta T_{\mu\nu}), \qquad (2.11)$$

where now the linearised Einstein operator takes the covariant form

$$\bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma} = \frac{1}{2}(\delta_{\mu}^{\rho} \delta_{\nu}^{\sigma} \bar{\nabla}^{2} - \delta_{\nu}^{\sigma} \bar{\nabla}_{\mu} \bar{\nabla}^{\rho} - \delta_{\mu}^{\sigma} \bar{\nabla}_{\nu} \bar{\nabla}^{\rho} + \bar{g}_{\mu\nu} \bar{\nabla}^{\sigma} \bar{\nabla}^{\rho} \\ + \bar{g}^{\rho\sigma} \bar{\nabla}_{\mu} \bar{\nabla}_{\nu} - \bar{s}_{\mu\nu} \bar{g}^{\rho\sigma} \bar{\nabla}^{2} - \bar{s}_{\mu\nu} \bar{R}^{\rho\sigma} + \delta_{\mu}^{\rho} \delta_{\nu}^{\sigma} \bar{R}), \qquad (2.12)$$

where $\bar{\nabla}_{\mu}$ and $\bar{R}$ are defined with respect to $\bar{s}_{\mu\nu}$. At this stage let us make an important remark on covariantisation of a theory whose form is known in flat space. Suppose we were given only the flat space Lagrangian (2.6) and asked to derive its covariant generalisation. Naïvely, one would replace all partial derivatives $\partial_{\mu}$ by covariant ones $\nabla_{\mu}$ and all occurrences of the Minkowski metric $\eta_{\mu\nu}$ by the more general background $\bar{s}_{\mu\nu}$. This procedure works for lower spin-particles coupled to gravity and results in well-known formulations of spin-0, spin-1/2 and spin-1 theories in curved space. For the spin-2 field, however, the procedure fails because naïve covariantisation of (2.6) does not result in the consistent form (2.12) obtained from linearising GR. In particular, ambiguities arise because the covariant derivatives do not commute and since there is no obvious guideline telling us which curvature terms to include in the linearised Einstein operator. On the other hand, knowledge of the nonlinear theory (2.7) allows us to straightforwardly arrive at the correct Lagrangian for a massless spin-2 field in curved background.

## 2.2. Massive spin-2 field in flat space

As has been known since the work of Fierz and Pauli from 1939 [3], the quadratic Lagrangian for a massive spin-2 excitation in flat space has the form
$$ \mathcal{L}_{\text{FP}} = \frac{1}{2} \left( h^{\mu\nu} \mathcal{E}_{\mu\nu}^{\rho\sigma} h_{\rho\sigma} - \frac{m_{\text{FP}}^2}{2} (h^{\mu\nu} h_{\mu\nu} - h^2) - \kappa h^{\mu\nu} T_{\mu\nu} \right), \quad (2.13) $$

where the kinetic operator is the same as in the massless case, given by (2.2), and we have included a source term $T_{\mu\nu}$. The equations of motion obtained from the Lagrangian for the massive particle are

$$ \mathcal{E}_{\mu\nu}^{\rho\sigma} h_{\rho\sigma} - \frac{m_{\text{FP}}^2}{2} (h_{\mu\nu} - \eta_{\mu\nu} h) = \frac{\kappa}{2} T_{\mu\nu}. \quad (2.14) $$

The mass term breaks the gauge invariance of the massless theory but, as a consequence of the linearised Bianchi identity, the divergence and trace of these equations give rise to five constraints

$$ \partial^{\mu} h_{\mu\nu} - \partial_{\nu} h = - \frac{\kappa}{m_{\text{FP}}^2} \partial^{\mu} T_{\mu\nu}, \quad (2.15a) $$

$$ h = \frac{\kappa}{3m_{\text{FP}}^2} T + \frac{2\kappa}{3m_{\text{FP}}^4} \partial^{\nu} \partial^{\mu} T_{\mu\nu}. \quad (2.15b) $$

The source is not necessarily conserved but, for simplicity, we shall anyway assume $\partial_{\mu} T^{\mu\nu} = 0$. This assumption certainly holds for any source that is derived from a diffeomorphism invariant matter coupling. For vanishing sources, the constraint equations imply that the massive spin-2 field is transverse and traceless.

Using (2.15) in the equations of motion, we can rewrite them as

$$ (\Box - m_{\text{FP}}^2) h_{\mu\nu} = \kappa \left( T_{\mu\nu} - \frac{1}{3} \left[ \eta_{\mu\nu} - \frac{1}{m_{\text{FP}}^2} \partial_{\mu} \partial_{\nu} T \right] \right) \quad (2.16) $$

and see that the massive spin-2 field satisfies a sourced Klein-Gordon equation. The constraints (2.15) remove five of the ten components in $h_{\mu\nu}$, leaving us with the five propagating degrees of freedom of a massive spin-2 particle.

An important observation made by Fierz and Pauli is that modifying the structure of the mass term, i.e. changing the numerical factor in front of $h^2$ in (2.13), introduces an additional degree of freedom into the theory. This happens because the trace constraint (2.15b) is lost and $h$ satisfies a dynamical equation of motion instead. It can furthermore be shown that the propagator corresponding to the extra degree of freedom comes with a residue of the wrong sign and therefore gives rise to a ghost instability<fnref n="7" /> [15]. This unwanted dynamical field is referred to as the *Boulware-Deser ghost*. It is this ghost mode that complicates the construction of a nonlinear interaction potential for a spin-2 field because, even when banned from the linear theory, the ghost notoriously reappears through the higher-order interactions.

On the other hand, tuning the linear mass potential to the Fierz-Pauli structure comes with its own problems that threaten the phenomenological viability of the theory as a modification of GR. Namely, as we already mentioned in the introduction, it was shown by vDVZ [9, 10] that the $m_{\text{FP}} \to 0$ limit of linear massive gravity does not continuously approach linearised GR. Nonlinear self-interactions for the spin-2 field may be able to cure this problem if they exhibit the Vainshtein mechanism [12]. Historically, this was one of the main motivations for the construction of a nonlinear theory of massive gravity. Another reason to search for a completion of the Fierz-Pauli mass term is the existence of the nonlinear closed form (2.7) in the massless theory.

<footnote n="7">⁷ See [79] for a detailed analysis of the propagator in Fierz-Pauli theory and its deformed version containing the ghost.</footnote>
Since consistency of linearised massive gravity in flat space requires tuning a coefficient in the mass potential, one can expect that a consistent (i.e. ghost-free) nonlinear potential cannot contain arbitrary interaction terms, but that the coefficients of certain terms will be related to each other by demanding the absence of the extra degree of freedom. If the fully nonlinear theory for massive spin-2 was known, it could be linearised around general backgrounds to give the covariantised version of (2.13). However, as in the massless case, it would be a very difficult task to derive the linear theory on arbitrary backgrounds by covariantising the flat-space Lagrangian. As we will see later in section 4.3, the most general linear theory for massive spin-2 can also be derived from a nonlinear action and has a rather complicated form.

### 3. Towards nonlinear spin-2 interactions

In the previous chapter we considered massless and massive theories for a linear fluctuation $h_{\mu\nu}$ in a fixed Minkowski background $\eta_{\mu\nu}$. In the nonlinear massless theory given by the Einstein-Hilbert action of GR, it was possible to combine the background and the fluctuation into a single nonlinear field $g_{\mu\nu}$. We would now like to do something similar in the massive case and construct a nonlinear self-interaction term without derivatives for the metric.

#### 3.1. General structure

Historically it was often assumed that any nonlinear theory for massive gravity must give rise to flat space solutions and a Fierz-Pauli mass term in the linear theory. On the other hand, it is well known that a consistent mass term may also be written down on maximally symmetric backgrounds<fnref n="8" /> [80–86], or even more generally on homogeneous and isotropic backgrounds [87–89]. This possibility motivates us to consider a general nonlinear theory without reference to flat Minkowski solutions and only demand that the correct number of degrees of freedom propagate at the nonlinear level (or, equivalently, around any background). Of course, when restricting to certain backgrounds the linearised version of the nonlinear theory must also reduce to the correct known structure.

What we call a mass term for a rank-2 tensor $g_{\mu\nu}$, must be a scalar density, i.e. it has to be a nontrivial scalar function $V(g)$ multiplied by the scalar density $\sqrt{g}$. Obviously, the scalar function $V(g)$ cannot have any loose covariant indices and, by definition, it should not contain any derivatives. But then, the only object at hand to contract the indices of the metric tensor $g_{\mu\nu}$ is the metric itself which, since $g^{\rho\mu}g_{\rho\nu} = \delta^{\rho}_{\nu}$, leads to a trivial cosmological constant contribution in the action. We conclude that there is no possibility to construct a covariant nonlinear interaction term for a spin-2 field using only one tensor field<fnref n="9" />. Hence, we are forced to introduce another field in order to build nonlinear contractions with $g_{\mu\nu}$. In principle, this could be any object with sufficient amount of indices, but the minimal choice is to work with a second rank-2 tensor which we shall call $f_{\mu\nu}$. The interaction potential will then be given by $\sqrt{g}$ multiplying a scalar function of $g^{\rho\mu}f_{\rho\nu}$. Note that, due to the existence of an additional tensor $f_{\mu\nu}$, we could in principle consider densitising by using for example $\sqrt{f}$ or $g^{1/4}f^{1/4}$. Since it is always possible to factor out $\sqrt{g}$ (for instance, by writing $\sqrt{f} = \sqrt{g} \det \sqrt{g^{-1}f}$ and

<footnote n="8">⁸ Interestingly, the vDVZ discontinuity turns out to be absent on de-Sitter backgrounds [80].</footnote>

<footnote n="9">⁹ Note the difference to e.g. the vector example, where the metric $g_{\mu\nu}$ can be used to contract the indices in nonlinear interaction terms for $A_\mu$ in a nontrivial way. Furthermore, in the massless spin-2 theory this problem did not occur because in the kinetic terms the indices of the metric could be contracted with derivative operators.</footnote>
absorbing the second factor into the potential) we may use that as the scalar density without any loss of generality.

In summary, we expect the nonlinear massive gravity action to be of the form

$$S_{\text{MG}} = m_g^2 \int d^4 x \sqrt{g} [R(g) - m^2 V(g^{-1}f)], \quad (3.1)$$

where $m$ is an arbitrary energy scale that sets the mass of the spin-2 field. In this setup, the second ‘metric’ $f_{\mu\nu}$ is a fixed background field that needs to be put into the theory by hand. From the viewpoint of field theory this is somewhat unusual because $f_{\mu\nu}$ is not determined by an equation of motion. In fact, there is no need to worry about this, since we will resolve this slightly confusing point later, when we introduce the fully dynamical bimetric theory that treats $g_{\mu\nu}$ and $f_{\mu\nu}$ on the same footing. For now, let us accept the possibility to work with a fixed reference metric $f_{\mu\nu}$ and investigate the consistency of this class of theories.

A simple example for a possible interaction term would be $V(g^{-1}f) = \text{Tr}(g^{-1}f)$. Interestingly, the corresponding action is closely related to another modification of gravity that goes under the name Eddington-inspired Born-Infeld theory [90]. Unfortunately, having only this term in the action gives rise to the Boulware-Deser ghost and is thus not a viable choice.

A necessary requirement on the interaction potential in (3.1) is that it reduces to the Fierz-Pauli mass term for parameter choices that permit a linearisation around flat space (i.e. when the reference metric $f_{\mu\nu}$ is flat and when flat background solutions for $g_{\mu\nu}$ exist), otherwise it will certainly propagate the Boulware-Deser ghost. But this requirement alone is not sufficient for consistency: a generic nonlinear interaction potential, even if it incorporates the Fierz-Pauli structure around flat space, will reintroduce the extra degree of freedom which leads to instabilities at the nonlinear level. This is anticipated from a simple degree-of-freedom counting in the full theory: the introduction of the interaction potential breaks the diffeomorphism invariance of GR. Therefore the four gauge symmetries are lost and the theory generically will propagate four additional degrees of freedom. Since $2+4=6$, this does not give the correct number for a massive spin-2 particle, but there is an extra degree of freedom in the theory. This is the nonlinear Boulware-Deser ghost and a constraint is needed to remove it from the spectrum of propagating modes. Tuning the linear mass term to the Fierz-Pauli combination ensures the presence of this constraint only in the linear theory.

We will see that, in order to obtain an additional constraint and thereby ensure the absence of the ghost beyond the linear level, it is necessary to impose strong restrictions on the structure of interactions. In fact, demanding the presence of a constraint will fix all but three coefficients of all possible interaction terms.

## 3.2. The Boulware-Deser ghost

In the construction of the consistent theory, we will follow the approach of Boulware and Deser who studied massive gravity in the Hamiltonian formulation [15, 91]. They claimed that, even though it is possible to remove the ghost from the spectrum of propagating modes at the linear level, it will return for any nonlinear interaction terms that are added to the Lagrangian (2.13). In order to show that this result is in fact incorrect, we first need to familiarise ourselves with variables suitable for the Hamiltonian formulation of GR.

### 3.2.1. ADM variables for GR.
The Hamiltonian formulation of GR traces back to the work by ADM from 1962 [60]. The authors decomposed the metric $g_{\mu\nu}$ into a scalar $N$ (lapse), a three-dimensional metric $\gamma_{ij}$, and a three-component vector $N_i$ (shift) as follows:
$$g_{\mu\nu} = \begin{pmatrix} -N^2 + N_i \gamma^{ij} N_j & N_j \\ N_i & \gamma_{ij} \end{pmatrix}, \qquad (3.2)$$

where $\gamma^{ij}$ denotes the inverse of $\gamma_{ij}$. This parametrisation essentially splits the metric into its time ($0\mu$) and spatial ($ij$) components. From (3.2) we can also compute the inverse of the metric

$$g^{\mu\nu} = \frac{1}{N^2} \begin{pmatrix} -1 & N^j \\ N^i & N^2 \gamma^{ij} - N^i N^j \end{pmatrix}. \qquad (3.3)$$

Here and in the following we raise the indices on the shift vector $N_i$ using the inverse spatial metric $\gamma^{ij}$.

It turns out that, in any theory with kinetic structure given by the Einstein-Hilbert term, the lapse $N$ and shift $N_i$ are non-dynamical gauge degrees of freedom because the Ricci scalar of $g_{\mu\nu}$ contains no derivatives on those fields. Hence, all of the propagating modes are contained in the spatial metric $\gamma_{ij}$ which has six independent components. The gauge invariance of GR further reduces the number of propagating degrees of freedom to two, along with the same number of corresponding canonical momenta.

More explicitly, in terms of the ADM variables the Einstein-Hilbert action (2.7) in vacuum becomes (up to a boundary term which, although crucial for certain applications, is not important for our considerations) [60]

$$S_{\text{GR}} = M_{\text{P}}^2 \int d^4 x (\pi^{ij} \partial_t \gamma_{ij} - N R^0 - N_i R^i), \qquad (3.4)$$

where, in terms of the curvature scalar $R_{(3)}$ of the metric $\gamma_{ij}$

$$R^0 = - \sqrt{\gamma} \left[ R_{(3)} + \gamma^{-1} \left( \frac{1}{2} (\pi^i)^2 - \pi^{ij} \pi_{ij} \right) \right], \\ R^i = -2 \partial_j \pi^{ij}. \qquad (3.5)$$

The conjugate momenta $\pi^{ij}$ of the six metric components $\gamma_{ij}$ are computed from the GR Lagrangian in the standard way which leads to expressions in terms of derivatives of the four-dimensional metric.

It is now evident that the action (3.4) does not contain dynamical terms for the scalar $N$ nor for the vector $N^i$. On top of that, these variables appear only linearly and therefore act as Lagrange multipliers whose equations of motion do not contain $N$ and $N^i$ themselves. This implies that these equations in fact correspond to four constraints $R^\mu = 0$, with $R^\mu = (R^0, R^i)$, on the remaining twelve variables $\gamma_{ij}$ and $\pi^{ij}$.

According to the theory of constrained Hamiltonian systems<fnref n="10" />, the presence of a gauge symmetry can now be seen in the Poisson algebra of the constraints, {$R^\mu, R^\nu$}. Here, the Poisson bracket for functions $A, B$ is defined as

$$\{A(x), B(y)\} \equiv \int d^3 z \left( \frac{\delta A(x)}{\delta \gamma_{ij}(z)} \frac{\delta B(y)}{\delta \pi^{ij}(z)} - \frac{\delta A(y)}{\delta \pi^{ij}(z)} \frac{\delta B(x)}{\delta \gamma_{ij}(z)} \right). \qquad (3.6)$$

<footnote n="10">¹⁰ For a review of this subject see e.g. [92].</footnote>
The result for the constraint algebra of GR reads [93, 94]

$$
\begin{align}
\{R^0(x), R^0(y)\} &= R^i(y) \frac{\partial}{\partial y^i} \delta^3(x-y) - R^i(x) \frac{\partial}{\partial x^i} \delta^3(x-y), \nonumber \\
\{R_i(x), R_j(y)\} &= R_i(y) \frac{\partial}{\partial y^j} \delta^3(x-y) - R_j(x) \frac{\partial}{\partial x^j} \delta^3(x-y), \nonumber \\
\{R^0(x), R_i(y)\} &= -R^0(y) \frac{\partial}{\partial x^i} \delta^3(x-y). \tag{3.7}
\end{align}
$$

Since these brackets are proportional to the constraints themselves, they vanish on the constraint surface, as they should in the presence of a gauge symmetry. In the language of Dirac, non-trivial gauge symmetries generate first class constraints. All four primary constraints will serve to put conditions on the remaining variables $\gamma_{ij}$ and $\pi^i$, while the Lagrange multipliers $N$ and $N^i$ remain undetermined at this stage. Since the Hamiltonian $H$ is itself a linear combination of $R$ and $R^i$, all constraints are automatically preserved in time since the time-evolution of any quantity is determined by $\frac{d}{dt}A(x) = \{A(x), H\}$.

After imposing the four conditions on $\gamma_{ij}$ and $\pi^i$, we can still use gauge transformations to remove another four degrees of freedom whose equations of motion will eventually determine $N$ and $N^i$. In total we therefore end up with four dynamical variables, corresponding to the two helicity states ($\pm 2$) of the massless graviton and their canonical momenta. This shows that, even nonlinearly, GR propagates the correct number of degrees of freedom for describing a massless spin-2 particle.

Note that throughout the above discussion we never actually wrote down the Hamiltonian, but remained in the Lagrangian formulation. In order to investigate the positivity of the energy<fnref n="11" /> this is naturally insufficient, but since here we were only interested in counting degrees of freedom there was no need to work directly with the Hamiltonian, which can trivially be obtained from (3.4).

### 3.2.2. The no-go theorem.
The ADM variables turn out to be very useful for investigating the consistency of massive gravity. We saw that at the linearised level the only consistent mass term is the one proposed by Fierz and Pauli given in (2.13). In order to demonstrate the presence of the ghost instability in the ADM formalism, in their paper [15] from 1972, Boulware and Deser studied the more general ‘mass term’ with arbitrary coefficient $a$

$$ \mathcal{L}_{\text{BD}} = \frac{1}{2} \left( h^{\mu\nu} \mathcal{E}_{\mu\nu}^{\rho\sigma} h_{\rho\sigma} - \frac{m^2}{2} (h^{\mu\nu} h_{\mu\nu} - ah^2) \right). \qquad (3.8) $$

The trace constraint $h = 0$ which removes the ghost mode does not exist for $a \neq 1$. Another way of seeing the additional mode is by noticing that, precisely for the Fierz-Pauli choice $a = 1$, the mass term is linear in the component $h_{00}$, which furthermore appears without time derivatives in the linearised Einstein operator. Thus, for $a = 1$, the equation of motion for $h_{00}$ is a constraint which removes one dynamical variable. For other values of $a$, the equation depends on $h_{00}$ itself in which case it does not constrain other components. As a direct consequence, a sixth degree of freedom is propagating.

We can make this more explicit using the variables introduced for the Hamiltonian analysis of GR in the previous subsection. From the ADM decomposition (3.2) we read off the decomposition of the fluctuation

<footnote n="11">¹¹ The positivity of the Hamiltonian has not been proven in the bimetric theory or in massive gravity in general, but has been studied for spherically symmetric configurations in [95].</footnote>
$$h_{\mu\nu} = g_{\mu\nu} - \eta_{\mu\nu} = \begin{pmatrix} 1 - N^2 + N_i \gamma^{ij} N_j & N_j \\ N_i & h_{ij} \end{pmatrix}, \qquad (3.9)$$

where $h_{ij} = \gamma_{ij} - \delta_{ij}$. Inserting this into the mass terms one finds

$$\frac{m^2}{4}(h^{\mu\nu}h_{\mu\nu} - ah^2) = -\frac{m^2}{4}(h^{ij}h_{ij} - a(h^i)^2 - 2N_iN_i + 2a(1 - N^2 + N_i\gamma^{ij}N_j)h^i_i + (1 - a)(1 - N^2 + N_i\gamma^{ij}N_j)^2). \qquad (3.10)$$

The last term gives rise to nonlinearities in $h_{00} = 1 - N^2 + N_i\gamma^{ij}N_j$ and vanishes only for $a = 1$. Moreover, when $h_{\mu\nu} = g_{\mu\nu} - \eta_{\mu\nu}$ is viewed as a small perturbation of $\eta_{\mu\nu}$, its ADM variables can be written as small fluctuations as well

$$\delta N = N - 1, \quad \delta N_i = N_i, \quad h_{ij} = \gamma_{ij} - \delta_{ij} \qquad (3.11)$$

and we can study their appearance to quadratic order in the mass term. For $a = 1$, the expression in (3.10) turns out to be linear in $\delta N$ at the quadratic level, such that the equation of motion of $\delta N$ gives rise to a constraint. On the other hand, we also see that the shift vector $\delta N_i$ does no longer appear only linearly which is of course consistent with the breaking of diffeomorphism symmetry by the mass term. Together with its associated secondary constraint, the $\delta N$ equation removes two degrees of freedom, one field plus its canonical momentum<fnref n="12" />. We end up with 12 - 2 = 10 degrees of freedom, describing the five polarisation states and corresponding conjugate momenta of the massive graviton.

Contrarily, for $a \neq 1$, there is a term involving $\delta N^2$ at the quadratic level. The constraint arising from the equation of motion of $\delta N$ is lost in that case because the equation now determines $\delta N$ itself instead of constraining other variables. There are thus 12 degrees of freedom, describing six propagating modes, one of which is the Boulware–Deser ghost.

The same situation occurs if we do not consider $h_{\mu\nu}$ as a small perturbation, i.e. look at more general backgrounds than $\eta_{\mu\nu}$. In that case, the ADM variables of $h_{\mu\nu}$ are no longer small fluctuations and we have to consider the theory beyond quadratic order. Boulware and Deser studied a class of corrections to the Fierz–Pauli mass term and found that those higher-order terms in $h_{\mu\nu}$ could never result in an expression that is linear in the lapse $N$. From this they concluded that the constraint that removes the ghost in the linear theory is destroyed and hence a theory of nonlinear spin-2 interactions can never be consistent [15, 91].

As pointed out in [57, 58] and as will be discussed in the next subsection, there exists a loophole in this no-go theorem for nonlinear massive gravity and it turns out that a consistent theory can be constructed<fnref n="13" />. In the following, we will deviate from the historical path and construct the consistent interaction potential directly in the redefined ADM variables that were used in the consistency proof of [59, 61] instead of presenting the derivation in the decoupling limit performed in [57, 58]. From our point of view this construction is the most efficient way to arrive at the action for nonlinear massive gravity and it has the further advantages that it (a) automatically ensures the absence of ghost in the full theory (i.e. away from the decoupling limit) and (b) immediately results in the generalised form of the action with arbitrary reference metric $f_{\mu\nu}$.

<footnote n="12">¹² The Fierz–Pauli theory is known to give rise to a secondary constraint coming from requiring the primary constraint to be preserved in time.</footnote>

<footnote n="13">¹³ Note also that the arguments given by Boulware and Deser are of perturbative nature. An idea that is rather different from everything we discuss here is that the negative norm states which plague the general higher-order interactions could be avoided non-perturbatively [96, 97].</footnote>
### 3.3. ADM variables for massive gravity

Our aim is to arrive at a nonlinear theory for massive spin-2 fields of the form (3.1) and, to this end, we shall discuss interactions in terms of ADM variables. For definiteness we will work in four-dimensions but all our considerations and conclusions generalise straightforwardly to any dimension. Before we start, let us briefly recapitulate the situation.

Since the kinetic term for the metric $g_{\mu\nu}$, in (3.1) is the same as in GR, the lapse and shift functions $N$ and $N^i$ will still appear without derivatives. However, the interaction potential will in general no longer be linear in these functions. Therefore, their equations of motion, instead of imposing constraints on the remaining variables, will now determine $N$ and $N^i$ themselves. The four gauge constraints are lost and, as explained in the previous subsection, the number of propagating modes will now generically be six, plus their corresponding canonical momenta. These are two phase-space degrees of freedom too many for the theory of a massive spin-2 field. Moreover, the Hamiltonian of a generic theory will not be positive definite [15, 91], signalling that the extra propagating mode is a ghost. A necessary requirement on any consistent interaction potential is therefore the presence of an additional constraint that removes the nonlinear Boulware–Deser ghost.

Before we set out to construct the potential featuring the constraint it should be noted that there are two other conditions which need to be fulfilled by a fully consistent theory. Firstly, the preservation of this constraint in time must itself provide a constraint in order to remove also the conjugate momentum and hence the full phase-space pair associated with the pathological degree of freedom. Secondly, the resulting Hamiltonian must be positive definite so that none of the surviving five spin-2 modes gives rise to an instability. It should be emphasised that positivity of the Hamiltonian of the nonlinear theories that we are about to discuss has never been proven in general and in fact seems not to be true without additional physical assumptions, see e.g. [95]. For instance, ghosts different from the Boulware–Deser mode may still propagate around certain backgrounds, an example being the Higuchi instability of the helicity-zero mode in de Sitter space [81]. In the literature, whenever massive gravity and bimetric theory are labelled as ‘consistent’ or ‘ghost-free’, one is usually referring only to the complete removal of the Boulware–Deser mode and its conjugate momentum. Throughout this review we frequently adhere to this conventional abuse of terminology.

#### 3.3.1. The loophole in Boulware and Deser's argument.
In order to investigate whether a particular structure in the potential $V(g^{-1}f)$ can give rise to a constraint and thus satisfy the first necessary condition on any consistent theory, we first decompose the second rank-2 tensor $f_{\mu\nu}$ into its own ADM variables<fnref n="14" />

$$f_{\mu\nu} = \begin{pmatrix} -L^2 + L_i \phi^{lk} L_k & L_j \\ L_i & \phi_{ij} \end{pmatrix}. \qquad (3.12)$$

Here, $\phi^{ij}$ denotes the inverse of the three-dimensional metric $\phi_{ij}$, $L_i$ is the shift-vector, and $L$ is the lapse of $f_{\mu\nu}$. We furthermore express the measure factor $\sqrt{g}$ in terms of the lapse $N$ and the determinant $\gamma$ of $\gamma_{ij}$

<footnote n="14">¹⁴ In fact, it is not automatically guaranteed that a simultaneous ADM split for $g_{\mu\nu}$ and $f_{\mu\nu}$ exists or, equivalently, that $N^2$ and $L^2$ are both positive definite. The assumption of simultaneous ADM decompositions, which shall be made here and in the following, is related to the existence of intersecting light cones for the two metrics. The details of this will be discussed in [98].</footnote>
$$ \sqrt{g} = N\sqrt{\bar{\gamma}}. \qquad (3.13) $$

With these expressions at hand, we can write a generic interaction potential in terms of ADM variables as<fnref n="15" />

$$ \sqrt{g} V(g^{-1}f) = N \sqrt{\bar{\gamma}} V(\gamma_{ij}, N, N^i; \phi_{ij}, L, L^i). \qquad (3.14) $$

According to the argument by Boulware and Deser, the right-hand side of this equation needs to be linear in the lapse $N$ in order to provide a constraint that removes the ghost mode. However, this is not entirely correct and this is exactly where the loophole lies: in fact, the necessary condition on the potential is slightly weaker because the constraint may be obtained after combining several of the equations for $N$ and $N^i$. In other words, if there exists a particular combination of these equations that is independent of $N^i$ and $N$, then this equation will constrain the remaining variables. At the level of the action, this means that we should allow for a field-dependent redefinition of the shift components, $N^i \to n^i$, that renders the Lagrangian linear in the lapse $N$. Let us explain this in a bit more detail. Suppose that a linear combination of the equations for $N$ and $N^i$ is independent of the lapse and thus corresponds to the constraint

$$ C = \frac{\delta S(N, N^j, \dots)}{\delta N} + C^i \frac{\delta S(N, N^j, \dots)}{\delta N^i} = 0, \qquad (3.15) $$

where the $C^i$ are some functions of the ADM variables. We can now make a redefinition of variables, $N^i(N, n^j, \dots) \equiv C^i_k(N, n^j, \dots)n^k$, such that the variation of the action with respect to the lapse becomes

$$ \begin{aligned} \frac{\delta S(N, n^j, \dots)}{\delta N} &= \frac{\delta S(N, n^j, \dots)}{\delta N} \Big|_{N^i} + \frac{\delta S(N, n^j, \dots)}{\delta N^i} \Big|_{N} \frac{\delta N^i(N, n^k, \dots)}{\delta N} \\ &= \frac{\delta S(N, n^j, \dots)}{\delta N} \Big|_{N^i} + \frac{\delta S(N, n^j, \dots)}{\delta N^i} \Big|_{N} \frac{\delta C^i_k(N, \dots)}{\delta N} n^k. \end{aligned} \qquad (3.16) $$

Here $|_{N^i}$ means that the function $N^i$ is kept fixed when the functional derivative is taken. From this we see that if we choose the redefined shift components $n^k$ such that they satisfy $\frac{\delta C^i_k(N, \dots)}{\delta N} n^k = C^i$ then the variation of the action with respect to the lapse gives precisely the constraint (3.15) which by assumption does not involve $N$. We shall thus look for a redefinition of the shift vector that renders the action linear in the lapse $N$.

Certainly, the redefinition, i.e. the matrix $C^i_k(N, \dots)$, must be linear in $N$ itself since the $N^i$ appear linearly in the kinetic term. Furthermore, since the redefined shift components $n^i$ are expected to appear in the constraint, they must be fully determined by their own equation of motion which therefore must not depend on $N$. To summarise, in order to fulfil the first necessary condition of obtaining a constraint we make two requirements on the potential in terms of the redefined shift $n^i$:

(i) linearity of the Lagrangian in $N$,

(ii) absence of $N$ in the $n^i$ equations of motion.

For now, we focus on the first requirement and we will see later that the potential that we construct by demanding only (i) automatically satisfies (ii).

### 3.3.2. Redefinition of ADM variables.
We now observe that the full interaction potential density, $\sqrt{g} V(g^{-1}f)$ in (3.1), already has a factor of $N$ in front coming from the measure

<footnote n="15">¹⁵ We will often make use of the notations $N^i \equiv \gamma^{ij}N_j$ and $L^i \equiv \phi^{ij}L_j$.</footnote>
factor (3.13). Hence, in order to satisfy requirement (i) listed above, the potential $V(g^{-1}f)$ written in redefined ADM variables must be of the form

$$V(g^{-1}f) = \frac{1}{N}V_1 + V_2, \qquad (3.17)$$

where $V_1$ and $V_2$ are functions only of $(\gamma_{ij}, n^i)$ (apart from the non-dynamical components of $f_{\mu\nu}$ which we choose to omit here), i.e. they are independent of $N$.

Recalling the ADM decomposition (3.2) of $g_{\mu\nu}$, we notice that $g_{\mu\nu}$ is quadratic in $N$ and in order to obtain inverse powers of $N$ we need to consider the inverse metric (3.3). The latter is quadratic in $1/N$ and the best we can achieve by a field redefinition which is linear in $N$ is to complete the dependence on $N$ into a perfect square such that taking a square-root can result in an expression linear in $1/N$. In other words, the only quantity that has a chance of giving something linear in $1/N$ after a linear redefinition of the $N^i$ is an object whose square is proportional to the inverse metric (3.3). We are thus led to consider a potential $V$ that is a function of the matrix $S \equiv \sqrt{g^{-1}f}$, defined via $S^2 = g^{-1}f$. This square-root matrix has a very nontrivial ADM decomposition and is certainly highly nonlinear in $1/N$ before any redefinition. However, we will now make use of the allowed redefinition of $N^i$ and demand that in terms of the new shift-vectors $n^i$ the square-root matrix $S$ is of the form [59, 61]

$$S = \frac{1}{N}A + B, \qquad (3.18)$$

where $A$ and $B$ are matrix-valued functions of $(\gamma_{ij}, n^i)$. The redefinition that leads to (3.18) as well as the explicit expressions for $A$ and $B$ can be obtained straightforwardly by the following method: square the right-hand side of the ansatz (3.18) and equate it with the ADM expression for $g^{-1}f$ obtained from (3.2) and (3.12), using the most general ansatz for the redefinition, $N^i = c^i + Nd^i$. Then compare the expressions on both sides order by order in $1/N$ to determine the vectors $c^i$ and $d^i$ in the shift redefinition as well as the matrices $A$ and $B$. This derivation was given in [61] and we discuss it in more detail in appendix B; here we simply state the result. The redefinition that renders the square-root matrix $S$ linear in $1/N$ takes a rather simple form [59, 61]

$$N^i = Ln^i + L^i + ND^i_k n^k. \qquad (3.19)$$

The $3 \times 3$ matrix $D$ on the right-hand side is a function of the variables $(\gamma_{ij}, n^i)$ and the non-dynamical spatial metric $\phi_{ij}$. Explicitly, it can be written in matrix notation as

$$D = \sqrt{\gamma^{-1}\phi} Q Q^{-1}, \qquad (3.20)$$

where we have defined another matrix $Q$ through

$$Q^i_j = x\delta^i_j + n^i n^k \phi_{kj}, \quad x = 1 - n^i \phi_{ij} n^j. \qquad (3.21)$$

Note that the definition (3.20) of the matrix $D$ that enters the redefinition involves a $3 \times 3$ square-root matrix. By introducing the shift vectors $n^i$ we have therefore reduced the dimension of the square-root matrix that appears in $S$ by one and, most importantly, simplified the dependence on the lapse $N$ which no longer appears under any square-root in (3.18). One may worry that a real solution for the $3 \times 3$ square-root in (3.20) does not always exist. However, we will show now that the variables can be further redefined to demonstrate the existence of real solutions for $D$ and, in fact, to remove the square-root matrix entirely.

### 3.3.3. On the existence of the redefinition.
The form of the redefinition (3.19) is not entirely unique. In fact, the original papers [59, 61] on massive gravity mainly worked with a set of
variables that slightly differs from the one presented here, whereas the choice of variables we made above is more suitable for application to bimetric theory [76]. Moreover, as was shown later in [99], it is possible to arrive at simpler expressions which are also more symmetric between the two metrics $g_{\mu\nu}$ and $f_{\mu\nu}$. In order to see this, let us decompose the two spatial metrics into ‘spatial dreibeins’

$$ \gamma_{ij} = e_i^a \delta_{ab} e_j^b, \qquad \phi_{ij} = \phi_i^a \delta_{ab} \phi_j^b. \qquad (3.22) $$

These expressions are invariant under rotations of the dreibeins, which means that $\tilde{\phi} = R\phi$ with $R^T = R^{-1}$ is an equivalent dreibein of $\phi$ and the rotations are a local symmetry of the theory. This freedom can be used to resolve the square-root appearing in $D$ (see (3.20)). Next, we redefine the shift vectors $n^i$ according to

$$ n^i = (\tilde{\phi}^{-1})^i_a \delta^{ab} v_b, \qquad (3.23) $$

where the new shifts $v^a$ carry a spatial Lorentz index and $\tilde{\phi} = R\phi$ is a gauge-fixed dreibein. The rotation matrix is chosen such that it satisfies

$$ \left[ \left( \hat{I} + \frac{1}{x + \sqrt{x}} v v^T \right) R \phi e^{-1} \right]^T = \left( \hat{I} + \frac{1}{x + \sqrt{x}} v v^T \right) R \phi e^{-1}, \qquad (3.24) $$

where $\hat{I}_{ab} = \delta_{ab}$. It is straightforward to obtain a solution for $R$ from this equation and this solution always exists, which follows directly from the polar decomposition theorem stating that any matrix can be symmetrised by a rotation. Inserting the new variables into (3.20) and using (3.24) to evaluate the square-root matrix, we find after some algebraic manipulations that $D$ takes the much simpler form

$$ D = \frac{1}{\sqrt{x}} e^{-1} \left( \hat{I} - \frac{1}{1 + \sqrt{x}} \hat{I}^{-1} v v^T \right) \tilde{\phi}, \qquad (3.25) $$

where we are using matrix notation and $\hat{I}_b^a = \delta_b^a$. Introducing the dreibeins and the new shift vectors has thus enabled us to get rid of any matrix square-root in the equations. The only square-root left is the scalar $\sqrt{x} = \sqrt{1 - v^T \hat{I}^{-1} v}$. This square-root has real solutions provided that the metric components satisfy the bound $v_a \delta^{ab} v_b < 1$. Interestingly, this bound has an interpretation connected with Lorentz transformations and it turns out that the redefined shift vector $v_a$ can be interpreted as a Lorentz velocity. We will come back to this point at the end of section 8.2.2. For now, let us assume that the bound is not violated and therefore the redefinition of shift vectors exists. In terms of the new variables, it takes the following simple and symmetric form [99]

$$ N^i = L^i + (L(\tilde{\phi}^{-1})^i_a + N(e^{-1})^i_a) \delta^{ab} v_b. \qquad (3.26) $$

Even though the structure of this expression is a little less complicated than (3.19), in order to stay closer to the conventions in the literature, we choose to continue using the original redefinition.

### 3.4. Construction of the ghost-free potential

As we show in appendix B, the matrices $A$ and $B$ in $\sqrt{g^{-1}f} = \frac{1}{N}A + B$ depend on the redefined ADM variables in the following way
$$ 
\begin{aligned}
\mathcal{A} &= \frac{1}{\sqrt{x}} \begin{pmatrix} L + n^k L_k & n^k \phi_{kj} \\ -(L + n^k L_k)(Ln^i + L^i) & -(Ln^i + L^i)n^k \phi_{kj} \end{pmatrix}, \\
\mathcal{B} &= \sqrt{x} \begin{pmatrix} 0 & 0 \\ D^i_j L^j & D^i_j \end{pmatrix},
\end{aligned} 
\quad (3.27) $$

in which $x = 1 - n^i \phi_{ij} n^j$ as before and the index on the shift vector $L_i$ is raised using the inverse spatial metric $\phi^{ij}$. Although the above expressions seem rather complicated, what is essential for the construction of the ghost-free potential is the structure of the matrices $\mathcal{A}$ and $\mathcal{B}$. Most importantly, $\mathcal{A}$ is a matrix of rank one and it can be written as the outer product of two vectors

$$ \mathcal{A} = uw^{\mathrm{T}}, \quad \text{where } u = \begin{pmatrix} 1 \\ -c^i \end{pmatrix} \quad \text{and} \quad w = \frac{1}{\sqrt{x}} \begin{pmatrix} a_0 \\ a_i \end{pmatrix}. \quad (3.28) $$

This will be the key property entering our constructive proof in the following. Having obtained the ADM expression with the desired dependence on $N$ for the square-root matrix $S = \sqrt{g^{-1}}f = \frac{1}{N}\mathcal{A} + \mathcal{B}$, we are ready to construct the interactions which give rise to an additional constraint that removes the Boulware–Deser ghost.

As explained above, the mass potential has to be a function of the square-root matrix $S = \sqrt{g^{-1}}f = \frac{1}{N}\mathcal{A} + \mathcal{B}$, multiplied by $\sqrt{g} = N\sqrt{\gamma}$. We assume this function to be, at least formally, expandable as a Taylor series in $S$. This commonly defines what one means by a matrix valued function anyway so is not really a serious restriction. Generically, this will give an expression that is nonlinear in the lapse $N$. The only way to ensure linearity in $N$ is to demand the absence of higher powers of $\frac{1}{N}\mathcal{A}$ in the expansion. Obviously, the simplest possible term is $\sqrt{g} \operatorname{Tr}(S) = \sqrt{\gamma} \operatorname{Tr}(\mathcal{A} + N\mathcal{B})$. At first sight it seems that all higher powers of $S = \frac{1}{N}\mathcal{A} + \mathcal{B}$ will involve higher powers of $1/N$. However, due to the special structure of the matrix $\mathcal{A}$ in (3.28) this is not quite correct and specific terms of higher order in $S$ can still be linear in $1/N$. Since $\mathcal{A}$ has rank one, it is a projection operator on a one-dimensional subspace. Owing to this property, there is a unique way of building polynomials of $\frac{1}{N}\mathcal{A} + \mathcal{B}$ that are linear in $\frac{1}{N}\mathcal{A}$. Namely, only an antisymmetric product of $\mathcal{A}$'s will automatically contain only one power of $\mathcal{A}$. Let us see how this works in detail by considering

$$ V(S) = \sum_{n=0}^{4} b_n \epsilon^{\mu_1 \mu_2 \dots \mu_n \lambda_{n+1} \dots \lambda_4} \epsilon_{\nu_1 \nu_2 \dots \nu_n \lambda_{n+1} \dots \lambda_4} S^{\nu_1}_{\mu_1} \dots S^{\nu_n}_{\mu_n}, \quad (3.29) $$

with arbitrary coefficients $b_n$ and totally antisymmetric tensors $\epsilon_{\mu\nu\rho\sigma}$. Since the simultaneous exchange of $\mu_i, \mu_j$ and $\nu_i, \nu_j$ only changes the sign of both $\epsilon$-tensors, it leaves the whole term invariant. Therefore, at $n$th order, the product of $S$'s under the sum can be written in the form

$$ \sum_{l=0}^{n} \binom{n}{l} \left(\frac{1}{N}\right)^l \mathcal{A}^{\nu_1}_{\mu_1} \dots \mathcal{A}^{\nu_l}_{\mu_l} \mathcal{B}^{\nu_{l+1}}_{\mu_{l+1}} \dots \mathcal{B}^{\nu_n}_{\mu_n}. \quad (3.30) $$

We insert this expression into (3.29) and obtain

$$ V(S) = \sum_{n=0}^{4} b_n \sum_{l=0}^{n} \binom{n}{l} \left(\frac{1}{N}\right)^l V_n(\mathcal{A}, \mathcal{B}), \quad (3.31) $$
with

$$
\begin{align}
V_n(\mathbf{A}, \mathbf{B}) &= \epsilon^{\mu_1\mu_2\cdots\mu_n} \lambda_{n+1}\cdots\lambda_4 \epsilon_{\nu_1\nu_2\cdots\nu_n} \lambda_{n+1}\cdots\lambda_4 \cdot \mathbf{A}^{\nu_1} \mu_1 \cdots \mathbf{A}^{\nu_l} \mu_l \mathbf{B}^{\nu_{l+1}} \mu_{l+1} \cdots \mathbf{B}^{\nu_n} \mu_n \nonumber \\
&= \epsilon^{\mu_1\mu_2\cdots\mu_n} \lambda_{n+1}\cdots\lambda_4 \epsilon_{\nu_1\nu_2\cdots\nu_n} \lambda_{n+1}\cdots\lambda_4 \cdot \mathbf{u}^{\nu_1} w_{\mu_1} \cdots \mathbf{u}^{\nu_l} w_{\mu_l} \mathbf{B}^{\nu_{l+1}} \mu_{l+1} \cdots \mathbf{B}^{\nu_n} \mu_n. \tag{3.32}
\end{align}
$$

Here we have used (3.28) in the second equality. Since all indices of the symmetric products of either $u^{\nu_i}$ or $w_{\mu_i}$ are contracted with the totally antisymmetric indices of the corresponding $\epsilon$-tensors, we find that in $V_n(\mathbf{A}, \mathbf{B})$ only terms with at most one $\mathbf{A}$ contribute. This implies that the sum over $l$ in (3.31) actually terminates at $l=1$ and hence $V(S)$ is linear in $1/N$, which is precisely the property that we were looking for. We therefore conclude that the most general form of the complete potential density which is linear in $N$ after the redefinition (3.19) is

$$
\sqrt{g} \sum_{n=0}^{4} b_n \epsilon^{\mu_1 \mu_2 \cdots \mu_n} \lambda_{n+1} \cdots \lambda_4 \epsilon_{\nu_1 \nu_2 \cdots \nu_n} \lambda_{n+1} \cdots \lambda_4 S^{\nu_1} \mu_1 \cdots S^{\nu_n} \mu_n \quad (3.33)
$$

and thus satisfies criterion (i) that we wrote down in section 3.3.1. In order to give rise to a
constraint it needs to also meet criterion (ii), that the $n^i$ equations following from the action
with the above potential need to be independent of the lapse $N$. As one can verify in a lengthy
but straightforward computation, this second requirement does not impose further restrictions
on the form of the potential but is automatically satisfied by (3.33). More explicitly, the $n^i$
equations are of the form [59, 61]

$$
\frac{\delta S}{\delta n^i} = \left( L \delta_i^j + N \frac{\delta (D_k^j n^k)}{\delta n^i} \right) E_j = 0, \quad (3.34)
$$

where $E_j$ does not involve $N$. Since the matrix multiplying $E_j$ is exactly equal to the Jacobian
$\frac{\delta N^j}{\delta n^i}$ of the redefinition, it is invertible by assumption (because otherwise the redefinition would
not be well-defined). We can thus multiply the equations by its inverse to arrive at the
equivalent equations $E_j = 0$ which do not depend on $N$.

We have thus derived the unique form of the interaction potential which gives rise to an additional constraint. Note in particular that there is only a finite number of terms giving a potential which is linear in N. Due to the antisymmetric structure of the interactions, the possible terms are limited by the dimension of spacetime. Before discussing their structure and properties in more detail in the next section, let us make some final remarks on the existence of the associated secondary constraint which arises from demanding the primary constraint to be preserved in time.

The secondary constraint: The above requirements that we used to construct the con-
sistent interaction potential were necessary but not sufficient: a secondary constraint is crucial
for the absence of the Boulware–Deser ghost because two constraints are needed to remove
both the ghost mode and its canonical momentum from the set of dynamical variables.

It was already motivated in [59, 61] that there is in fact a secondary constraint arising
from demanding the primary constraint to be constant in time. In order to compute the time
evolution of the constraint one uses the Poisson bracket (3.6) and the Hamiltonian, which for
the massive gravity action (3.1) with potential (3.33) is of the form

$$
H = \int d^3 y (\mathcal{H}_0 - N C). \tag{3.35}
$$

Here $\mathcal{H}_0$ is independent of $N$ and $C$ is the Hamiltonian constraint obtained from varying the
action with respect to $N$. The time evolution of $C$ then reads
$$ \frac{d\mathcal{C}(x)}{dt} = \{\mathcal{C}(x), H\} = 0. \qquad (3.36) $$

As in the case of the primary constraint, this equation should be independent of N, because otherwise it would determine N instead of constraining $\gamma_{ij}$ and $\pi^{ij}$. Inserting (3.36) into (3.35) gives

$$ \int d^3x \left( \{\mathcal{C}(x), \mathcal{H}_0(y)\} - N \{\mathcal{C}(x), \mathcal{C}(y)\} \right) = 0. \qquad (3.37) $$

Since $\mathcal{C}$ and $\mathcal{H}_0$ are independent of N, we need the Poisson bracket among the constraints $\{\mathcal{C}(x), \mathcal{C}(y)\}$ to vanish. Before this was actually demonstrated to be the case, there had been objections against the theory claiming that $\{\mathcal{C}(x), \mathcal{C}(y)\}$ could be nonzero [100]. The issue was resolved when the secondary constraint was eventually shown to exist in a detailed calculation in [62]. We will not repeat this analysis here but instead refer the interested reader to the original reference as well as other subsequent independent confirmations of these results [27, 63, 67–69].

# 4. Ghost-free nonlinear massive gravity

In the previous section we derived the nonlinear interaction potential for massive gravity equipped with an additional constraint that removes the Boulware–Deser ghost. Here we will study its properties in more detail and discuss the subclass of dRGT models as well as the gauge invariant Stückelberg formulation.

## 4.1. Most general action

For notational purposes, a more convenient way of writing the consistent potential in (3.33) is to introduce the elementary symmetric polynomials $e_n(S)$ of the matrix $S$. These can be defined in terms of totally antisymmetric tensors (with unit weight)

$$ e_n(S) = \frac{1}{n!(4-n)!} \epsilon^{\mu_1\mu_2\cdots\mu_n} \lambda_{n+1}\cdots\lambda_4 \epsilon_{\nu_1\nu_2\cdots\nu_n} \lambda_{n+1}\cdots\lambda_4 S^{\nu_1}_{\mu_1}\cdots S^{\nu_n}_{\mu_n}, \qquad (4.1) $$

with $e_0(S) \equiv 1$. More properties of the elementary symmetric polynomials as well as precise definitions of the $\epsilon$-tensors as anti-symmetrisation operators are summarised in appendix A. In terms of these the complete action for ghost-free massive gravity with general reference metric takes the form

$$ S_{\text{MG}} = m_g^2 \int d^4x \sqrt{g} \left[ R(g) - 2m^2 \sum_{n=0}^{4} \beta_n e_n(S) \right], \qquad (4.2) $$

with $S = \sqrt{g^{-1}}f$, Planck mass $m_g$ and spin-2 mass scale $m$. Furthermore, we have introduced the rescaled coefficients $\beta_n = b_n n!(4-n)!/2$ for $n = 0, \dots, 4$. Out of these five parameters, only three are truly measuring interaction strengths: since $e_0(S) = 1$, the $\beta_0$-term is simply a cosmological constant for the dynamical metric $g_{\mu\nu}$. Moreover, the last term in the sum that is proportional to $\beta_4$ is just a cosmological constant term for $f_{\mu\nu}$ because $e_4(S) = \det S$ and hence $\sqrt{g} e_4(S) = \sqrt{f}$. This term is therefore independent of $g_{\mu\nu}$ and does not contribute to the equations of motion. Nevertheless, we choose to include it in the action because it will become relevant when we give dynamics to $f_{\mu\nu}$ later on.

The above action is a nontrivial generalisation of the dRGT model [57, 58] which we shall discuss in the next subsection. Its above form (with general reference metric $f_{\mu\nu}$, finite sum over n and in terms of the elementary symmetric polynomials) was first presented in [46].
The equations of motion for $g_{\mu\nu}$ obtained from (4.2) are

$$R_{\mu\nu}(g) - \frac{1}{2}g_{\mu\nu}R(g) + m^2 V_{\mu\nu}^g(g, f; \beta_n) = 0, \quad (4.3)$$

where the first two terms are the usual Einstein tensor while the contribution from the interaction potential is

$$V_{\mu\nu}^g = g_{\mu\rho} \sum_{n=0}^{3} (-1)^n \beta_n (Y_n)^\rho_\nu(S), \quad (4.4)$$

where we have defined the matrix functions

$$(Y_n)^\rho_\nu(S) \equiv \sum_{k=0}^{n} (-1)^k e_k(S) (S^{n-k})^\rho_\nu. \quad (4.5)$$

There is no equation of motion for $f_{\mu\nu}$ whose form therefore needs to be put in by hand<fnref n="16" />. Taking the divergence of the above equations and using the Bianchi identity $\nabla^\mu G_{\mu\nu} = 0$ satisfied by the Einstein tensor $G_{\mu\nu}$, we arrive at a set of Bianchi constraints

$$\nabla^\mu V_{\mu\nu}^g = 0. \quad (4.6)$$

These remove four degrees of freedom while the remaining extra scalar (the Boulware–Deser ghost) is eliminated by the additional constraint present in the special structure of (4.2). A covariant expression for this scalar constraint is difficult, if not impossible, to obtain in general. Its explicit form for certain regions in the parameter space is provided in [101, 102] which make use of the vierbein formulation that we shall discuss in section 8.2. For the same restricted parameter choices, it is also possible to identify a covariant constraint in the linear theory around arbitrary backgrounds, see section 4.3.

## 4.2. dRGT theory

The theory first derived by dRGT in [57, 58] was defined for flat reference metric $f_{\mu\nu} = \eta_{\mu\nu}$. The original construction uses Stückelberg fields (see section 4.4) and its ‘ghost proof’ is valid only in the scalar sector of a decoupling limit that strongly relies on the flat reference metric<fnref n="17" />. In contrast, the construction that we have presented here is based on the results of [59, 61] and is valid for all reference metrics. It demonstrates the consistency of the full nonlinear theory away from any limiting approximation.

It is worth discussing the dRGT model and its relation to the Hassan–Rosen formulation in a bit more detail. Using the observation that the sum in the interaction potential terminates and can be given in terms of elementary symmetric polynomials [46], the dRGT action can be written in the form [57, 58] (using the conventions of [2])

$$S_{\text{dRGT}} = m_g^2 \int d^4 x \sqrt{g} \left( R(g) + \frac{m^2}{2} \sum_{n=0}^{4} n! (4-n)! \alpha_n e_n(\mathcal{K}) \right), \quad (4.7)$$

<footnote n="16">¹⁶ Imposing an equation of motion for the reference metric, one can integrate out $f_{\mu\nu}$ from the action but in this case the theory becomes dynamically equivalent to GR.</footnote>

<footnote n="17">¹⁷ For the sake of historical accuracy, let us emphasise again, that [57, 58] also pointed out the loophole in Boulware and Deser’s argument. Moreover, the lowest-order terms of the shift redefinition in a perturbative expansion around flat space were computed, but the results do not agree with the nonlinear result (3.19) derived in [59, 61].</footnote>
where

$$ K_{\nu}^{\mu} = \delta_{\nu}^{\mu} - (\sqrt{g^{-1}\eta})_{\nu}^{\mu}, \quad \alpha_{0} = \alpha_{1} = 0, \quad \alpha_{2} = 1. \qquad (4.8) $$

We can arrive at this action starting from (4.2) by setting $f_{\mu\nu} = \eta_{\mu\nu}$ and taking

$$ \beta_n = - \frac{(-1)^n}{4} \sum_{k=n}^{4} \binom{4-n}{k-n} k!(4-k)! \alpha_k. \qquad (4.9) $$

This follows from the identity (A.4) satisfied by the elementary symmetric polynomials. Hence the consistency of massive gravity with general $f_{\mu\nu}$ implies the absence of ghost in the dRGT action. In contrast, there is no obvious way of getting to (4.2) from (4.7) and results obtained in the latter do not generalise automatically to the former.

Let us briefly explain why the lowest-order $\alpha_n$ parameters in the dRGT action are fixed while they remain arbitrary in the Hassan–Rosen formulation, since the exact reason for this is sometimes obscured. For $g_{\mu\nu} = \eta_{\mu\nu}$ to be a solution to the equations of motion following from (4.7), the cosmological constant for $g_{\mu\nu}$ must vanish. It is straightforward to verify that this<fnref n="18" /> requires $4\alpha_0 + \alpha_1 = 0$, since this combination is proportional to the effective cosmological constant for $g_{\mu\nu}$. Next, in order to remove terms linear in the perturbation $h_{\mu\nu} = g_{\mu\nu} - \eta_{\mu\nu}$ in the quadratic action (so-called ‘tadpoles’), one must<fnref n="19" /> enforce $\alpha_1 = 0$. Finally, the interaction parameters contain a redundant overall scale that can be absorbed into $m^2$. One way to get rid of this redundancy is to demand that $m^2$ corresponds to the squared Fierz–Pauli mass in the quadratic action around flat space. This requirement finally fixes $\alpha_2 = 1$.

The above action (4.7) thus gives rise to Minkowski solutions for $g_{\mu\nu}$ and linearising the theory around these backgrounds gives precisely the Fierz–Pauli Lagrangian (2.13) with $m_{\text{FP}}^2 = m^2$. In this sense the dRGT theory can be viewed as the consistent nonlinear completion of the Lagrangian for a massive spin-2 field in flat space. In the following, we will turn to the linear theory around arbitrary background solutions of the more general theory (4.2).

### 4.3. Linear theory on arbitrary background

As we pointed out in section 2.2, naïvely covariantising the linear theory for a massive spin-2 field in flat background does not result in a consistent action. Attempts to find the correct equations describing five helicity states around any background were made already before the ghost-free nonlinear theory was known. For instance, the authors of [103] were able to write down the consistent linearised theory to first order in a small-curvature expansion.

On the other hand, the knowledge of the full nonlinear action (4.2) that avoids the ghost mode makes it possible to derive the linear theory around any background solution<fnref n="20" />. This computation was first carried out in [104, 105] (see also [106–108]). The analysis is complicated by the presence of the square-root matrix $\sqrt{g^{-1}f}$ whose perturbation around general

<footnote n="18">¹⁸ It should be noted that this first requirement can actually be avoided by demanding only that $g_{\mu\nu} = c^2\eta_{\mu\nu}$ is a solution. This is still flat of course and fixes c instead of a parameter of the action.</footnote>

<footnote n="19">¹⁹ From the bimetric perspective, this second requirement would mean that the effective cosmological constant for $f_{\mu\nu}$ vanishes since the latter is proportional to $\alpha_1 \sim \beta_1 + 3\beta_2 + 3\beta_3 + \beta_4$, see equation (5.11). The first requirement of vanishing cosmological constant for $g_{\mu\nu}$ can also be written as $4\alpha_0 + \alpha_1 \sim \beta_0 + 3\beta_1 + 3\beta_2 + \beta_3 = 0$.</footnote>

<footnote n="20">²⁰ Note that in order to work with the most general backgrounds, we have to start from nonlinear massive gravity in the generalised formulation (4.2), where $f_{\mu\nu}$ is entirely arbitrary. The set of solutions for the metric $g_{\mu\nu}$ in the dRGT model (4.7) with flat reference metric is much smaller and hence does not allow us to derive the most general form for the linearised theory. For instance, it is not possible to obtain the quadratic action for massive spin-2 in de Sitter space from (4.7).</footnote>
backgrounds is rather nontrivial due to the matrices $g_{\mu\nu}$ and $f_{\mu\nu}$ in general being non-commuting. The full expression can however be obtained using the Cayley–Hamilton theorem for matrices. An alternative way of arriving at the linearised theory is to redefine the dynamical fluctuation variables in order to simplify the square-root variation [108].

The general structure of the quadratic Lagrangian for perturbations $\delta g_{\mu\nu}$ around arbitrary backgrounds $\bar{g}_{\mu\nu}$ is

$$ \mathcal{L} = \delta g_{\mu\nu} \mathcal{E}^{\mu\nu\rho\sigma} \delta g_{\rho\sigma} + \frac{m^4}{m_g^2} \delta g_{\mu\nu} \mathcal{V}^{\mu\nu\rho\sigma} \delta g_{\rho\sigma}, \quad (4.10) $$

where $\mathcal{E}^{\mu\nu\rho\sigma}$ is a function of the background metric $\bar{g}_{\mu\nu}$ and quadratic in its associated covariant derivative $\bar{\nabla}^{\mu}$, whereas $\mathcal{V}^{\mu\nu\rho\sigma}$ depends on the background metric, its curvatures and the reference metric $f_{\mu\nu}$. The corresponding equations of motion for the fluctuations read

$$ \delta E_{\mu\nu} \equiv \mathcal{E}_{\mu\nu}^{\rho\sigma} \delta g_{\rho\sigma} + \frac{m^4}{m_g^2} \mathcal{V}_{\mu\nu}^{\rho\sigma} \delta g_{\rho\sigma} = 0. \quad (4.11) $$

Using the background equations (4.3), it is always possible to convert curvatures of $\bar{g}_{\mu\nu}$ into functions of the background and reference metric. Eliminating the reference metric through the background equations is more difficult since the equations are nonlinear in $f_{\mu\nu}$. In general it is not possible to obtain a closed form for the solution. However, quite remarkably, it turns out that if $\beta_2 = \beta_3 = 0$ a closed form solution exists and the equations can be used to get rid of all appearances of the reference metric. Then (4.11) describes five massive spin-2 degrees of freedom propagating in an arbitrary background given solely by $\bar{g}_{\mu\nu}$ and its curvature $\bar{R}_{\mu\nu}$. In the more general case it is still possible to solve the background equations perturbatively in curvatures and reproduce the results of [103].

The absence of the Boulware–Deser ghost in the linearised equations follows from the existence of the additional constraint in the nonlinear theory from which (4.10) was obtained. In order to bypass the ADM analysis and show explicitly that the above equations contain the same number of propagating degrees of freedom as in the Fierz–Pauli case (2.14) in flat space, one needs to identify the analogues of the constraints (2.15). The vector constraints, i.e. the generalisation of $\partial^{\mu}h_{\mu\nu} - \partial_{\nu}h = 0$, are easily obtained by taking the (covariant) divergence of the equations and using the linearised version of the Bianchi identity. The additional scalar constraint which removes the ghost is more difficult to identify. For models with $\beta_3 = 0$, it has been found in [104, 105, 108]. It corresponds to the following combination of equations

$$ \frac{1}{m^2} (S^{-1})^{\nu}_{\rho} \nabla^{\rho} \nabla^{\mu} \delta E_{\mu\nu} + \frac{\beta_1}{2} \bar{g}^{\mu\nu} \delta E_{\mu\nu} + \beta_2 S^{\nu}_{\rho} \bar{g}^{\rho\mu} \delta E_{\mu\nu} = 0, \quad (4.12) $$

in which all terms containing two derivatives on $\delta g_{\mu\nu}$ cancel out. This equation is the analogue of $h = 0$, which it reduces to when the background is an Einstein spacetime. When $\beta_3 \neq 0$, it seems to be impossible to obtain the constraint in a covariant way [102, 108] and the reason for that is not yet fully understood.

### 4.4. Gauge invariant massive gravity

A common perspective is to take the fixed reference metric $f_{\mu\nu}$ to be fully specified in a given coordinate frame which implies that the invariance under GCTs of GR is broken by the mass term for the graviton. This is because the metric $g_{\mu\nu}$ itself transforms as a covariant rank-two tensor under GCTs while the fixed background metric $f_{\mu\nu}$ does not if it is taken to be fixed in a given coordinate system. Therefore, in this view, objects like $g^{\rho\mu}f_{\mu\nu}$ do not transform as
tensors under diffeomorphisms. Of course covariance is not a true symmetry and it is always possible to simply view $f_{\mu\nu}$ as a fixed reference geometry but still allow it to transform as a tensor under GCTs, e.g. to be flat but not necessarily on the cartesian Minkowski form. A constructive way of restoring any gauge symmetry in the action is the so called Stückelberg trick: by introducing new gauge degrees of freedom, so-called Stückelberg fields, that transform in a certain way under the gauge group, one can rewrite the theory in a manifestly gauge invariant way while keeping track of the thereby introduced redundant degrees of freedom. In turn, fixing the so-called physical (or unitary) gauge will give back the original action.

**4.4.1. The stückelberg trick.** One possibility to reintroduce diffeomorphism invariance in massive gravity is to perform a gauge transformation on the dynamical metric $g_{\mu\nu}$ and treat the gauge parameters as new fields in the action. Under a GCT, $x^\alpha \mapsto x'^\alpha(x)$, the metric transforms as

$$g_{\mu\nu}(x) \mapsto \frac{\partial x'^\alpha}{\partial x^\mu} \frac{\partial x'^\beta}{\partial x^\nu} g_{\alpha\beta}(x'). \qquad (4.13)$$

After performing this GCT, one interprets the $x'^\alpha$ as dynamical fields. Another equivalent way, that turns out to be more convenient for studying massive gravity, is to mimic this kind of transformation on the reference metric $f_{\mu\nu}$, which is then taken to not transform under GCTs. But since a simultaneous coordinate transformation of $g_{\mu\nu}$ and $f_{\mu\nu}$ would be a symmetry of the massive gravity action, instead of transforming $g_{\mu\nu}$ we can choose to perform the transformation on $f_{\mu\nu}$. That is to say, we replace $f_{\mu\nu}$ by its covariantised form

$$f_{\mu\nu} \rightarrow \frac{\partial \varphi^\alpha}{\partial x^\mu} \frac{\partial \varphi^\beta}{\partial x^\nu} \bar{f}_{\alpha\beta}. \qquad (4.14)$$

Here, $\bar{f}_{\alpha\beta}$ is a new fixed background metric that does not transform under the gauge group. The ‘coordinates’ $\varphi^\alpha$ are promoted to dynamical fields that transform as scalars under GCTs. Introducing these Stückelberg fields ensures that objects such as

$$\mathrm{Tr}(g^{-1}f) = g^{\mu\nu} \frac{\partial \varphi^\alpha}{\partial x^\mu} \frac{\partial \varphi^\beta}{\partial x^\nu} \bar{f}_{\alpha\beta} \qquad (4.15)$$

now are manifestly diffeomorphism invariant. In addition to the equations of motion for the metric $g_{\mu\nu}$, we now also consider the $\varphi^\alpha$ equations and, as we will see below, this system of equations is equivalent to the original $g_{\mu\nu}$ equation without Stückelberg fields. This method of restoring gauge invariance, in contrast to treating $f_{\mu\nu}$ as a fixed geometry tensor, is particularly useful for considering so called decoupling limits of a theory: limits/regimes where certain operators become dominant and physical modes decouple.

**4.4.2. The decoupling limit.** Let us take $\bar{f}_{\alpha\beta} = \eta_{\alpha\beta}$ and expand the Stückelberg fields around the identity transformation, which corresponds to considering infinitesimal GCTs

$$\varphi^\alpha = x^\alpha - \pi^\alpha. \qquad (4.16)$$

Furthermore, let us decompose the perturbations into a transverse and a longitudinal mode

$$\pi^\alpha = \frac{\hat{\pi}^\alpha}{mm_g} + \frac{\partial^\alpha \pi}{m^2 m_g}, \qquad (4.17)$$
with $\partial_{\alpha}\hat{\pi}^{\alpha} = 0$. The ghost-free dRGT potential (4.7) was originally constructed using the gauge invariant formulation and performing a scaling limit to separate the interactions for the longitudinal component $\pi$.

The power of the Stückelberg formalism lies in a conjecture made in [47], stating that the ghost instability of nonlinear massive gravity can be traced back to higher-derivative interactions of the $\pi$ fields in flat space, where $\tilde{f}_{\alpha\beta} = \eta_{\alpha\beta}$. This assumption relies on an analogy to the Goldstone equivalence theorem for spin-1 fields [49, 109]. In that case, the longitudinal modes of the gauge fields have been shown to carry all information needed for computing scattering processes at high energies. Although the theorem has not been proven for the spin-2 case, it is reasonable to start with the assumption of its validity in order to investigate the stability of massive gravity in a certain parameter limit defined through

$$m_g \rightarrow \infty, \quad m \rightarrow 0, \quad \Lambda_3 \equiv (m^2 m_g)^{\frac{1}{3}} = \text{const.} \qquad (4.18)$$

The metric is decomposed into flat background $\eta_{\mu\nu}$ and fluctuations $h_{\mu\nu}$ according to $g_{\mu\nu} = \eta_{\mu\nu} + \frac{h_{\mu\nu}}{m_g}$. The fields in the Stückelberg decomposition (4.17) already have the correct normalisations and, in a first approach, the vector modes are set to zero. Then the decomposition of the tensor $g^{-1}f$ that appears in the interaction potential into flat background and perturbation reads, in matrix notation

$$g^{-1}f = \left(1 + \frac{\eta^{-1}h}{m_g}\right)^{-1} \left(1 + \frac{\Pi}{\Lambda_3^{\frac{2}{3}}}\right)^2, \quad \Pi^{\mu}_{\nu} \equiv \eta^{\mu\rho} \partial_{\rho} \partial_{\nu} \pi. \qquad (4.19)$$

The conjecture of [47] suggests now to construct the interaction potential in such a way that there appear no higher-order interactions of the longitudinal modes. Since each of these modes always comes with two derivatives, their interactions generically give rise to higher-derivative terms These can be shown to lead to ghost instabilities due to the famous Ostrogradsky theorem [110] and must therefore be avoided in a consistent theory (see [111] for a recent discussion of this theorem). dRGT managed to tune the coefficients of the $\partial_{\mu}\partial_{\nu}\pi$ interaction terms in such a way that the higher-derivative terms combined into total derivatives that could be dropped from the action [57, 58]. Interestingly, the resulting action for the longitudinal modes (after the kinetic terms have been diagonalised) resembles precisely the Galileon interactions [112]. In fact, from (4.19) it is a trivial observation to see that by taking the square-root the expression becomes linear in $\Pi_{\mu\nu}$ and then any antisymmetric products of this square-root will at most carry two derivatives on any single $\pi$.

A few more comments are in order. It is sometimes argued (see e.g. [2, 113, 114]) that any massive gravity theory can be described as a theory of GR plus four scalar fields. In other words, this would imply that it is always possible to reduce the number of degrees of freedom of the massive fluctuation $h_{\mu\nu}$ to two (corresponding to a massless spin-2 mode) and let the remaining three reside in the Stückelberg fields. Let us briefly clarify why this picture is not quite correct. Firstly, an equivalent statement would be that it is always possible to render the background solution for $g_{\mu\nu}$ flat by performing a coordinate transformation. As already pointed out in [115], this is obviously not correct since the flat space action and equations are not equivalent to the ones in curved space. Secondly, the counting of degrees of freedom in GR is an on-shell statement and only holds for the massless equations of motion, e.g. $\Box h_{\mu\nu} = 0$ on a flat background, which do not arise in a massive gravity theory. Therefore, the number of degrees of freedom in $h_{\mu\nu}$ cannot be reduced to two by simply performing a gauge transformation. Thirdly, from the transformation properties of the $\pi^{\alpha}$ fields in (4.16) under linearised coordinate transformations ($\delta_{\xi}\pi^{\alpha} = -\xi^{\alpha}$), it is clear that these fields can only mix with the four gauge modes of $g_{\mu\nu}$. In other words, they are gauge trivial by
construction. Physical degrees of freedom can therefore not be fully transferred to the gauge modes because otherwise they would not contribute to interactions between conserved sources (for which $\nabla_{\mu} T^{\mu\nu} = 0$).

Although these caveats cast some serious doubts on the conclusiveness of the consistency analysis in the decoupling limit (other objections were raised in e.g. [116]), its historical importance for the discovery of ghost-free massive gravity should of course not be underestimated. The complete expression for the dRGT theory in the decoupling limit, including the transverse vector modes in (4.17), was later worked out in [117]. The limit is furthermore useful for understanding the validity regime of the effective field theory and the identification of the strong-coupling scale of massive gravity. We refer the interested reader to the review [2] where these subjects are discussed in more detail.

### 4.4.3. Absence of ghost in the Stückelberg formulation.
While the proofs of absence of the Boulware–Deser ghost given in [59, 61] were formulated in the gauge fixed version of massive gravity, i.e. without introducing the Stückelberg fields, it is expected that restoring the gauge invariance does not alter the dynamics of the theory. In fact, since gauge symmetry is merely a redundancy in the description of the underlying physics, the two formulations are completely equivalent. In particular, absence of ghost in the gauge fixed version of massive gravity implies that the Stückelberg formulation cannot exhibit the instability either. Nevertheless, in [118] it was claimed that the ghost generically reappears in the Stückelberg sector and the author of [119, 120] could not find a constraint to remove the ghost. First attempts to disprove these statements were made in [114, 121], before the absence of ghost was conclusively shown in [63].

The fact that the dynamics remain unaltered in the gauge invariant theory can be demonstrated by considering the equations of motions for the Stückelberg fields. After making the replacement (4.14) for $f_{\mu\nu}$, the mass potential $V(g, f)$ becomes a function of the scalar fields and their equations of motion read

$$ \frac{\delta}{\delta\varphi^{\alpha}} V(g, \varphi) = 0. \qquad (4.20) $$

However, this equation will not give rise to any new dynamics because it is, in fact, already implied by the Bianchi constraint

$$ \nabla^{\mu} V_{\mu\nu}(g, \varphi) = 0, \qquad (4.21) $$

where

$$ V_{\mu\nu}(g, \varphi) \equiv \frac{-2}{\sqrt{g}} \frac{\delta}{\delta g^{\mu\nu}} (\sqrt{g} V(g, \varphi)). \qquad (4.22) $$

In order to see the equivalence, consider a gauge transformation $\delta x^{\mu} = \xi^{\mu}$ of the action (4.2) involving Stückelberg fields

$$ \delta S = m_g^2 \int d^4 x \sqrt{g} \left[ (G_{\mu\nu} + m^2 V_{\mu\nu}) \delta g^{\mu\nu} - 2m^2 \frac{\delta V}{\delta \varphi^\alpha} \delta \varphi^\alpha \right], \qquad (4.23) $$

where $G_{\mu\nu}$ is the Einstein tensor and the variations of the fields under the gauge transformation are $\delta g^{\mu\nu} = 2\nabla^{(\mu}\xi^{\nu)}$ and $\delta \varphi^\alpha = -\xi^\mu \partial_\mu \varphi^\alpha$. The action is now invariant under GCTs, $\delta S = 0$. Thus integrating by parts and using the Bianchi identity, $\nabla^\mu G_{\mu\nu} = 0$, we can derive the following identity
$$ \nabla^{\mu} V_{\mu\nu} = \frac{\delta V}{\delta \varphi^{\alpha}} \partial_{\nu} \varphi^{\alpha}. \qquad (4.24) $$

Since the Stückelberg fields were introduced to mimic nonsingular coordinate transformations, the matrix $\partial_{\nu}\varphi^A$ is invertible and hence the Bianchi constraint $\nabla^{\mu}V_{\mu\nu} = 0$ is equivalent to the $\varphi^\alpha$ equations of motion (4.20). As expected, the redundant gauge degrees of freedom $\varphi^\alpha$ therefore do not introduce new dynamics into the theory.

From these general arguments it is obvious that the consistency proof of [59, 61] has to be valid for both the gauge fixed and the gauge invariant version of massive gravity. This was eventually confirmed in a Hamiltonian analysis where the constraint was shown to exist in the massive gravity action (with $\beta_2 = \beta_3 = 0$) including Stückelberg fields [63].

### 4.5. Potential shortcomings of massive gravity

In order to conclude the presentation of nonlinear massive gravity, we list a few drawbacks of the theory in its present form. These can be viewed as motivations to look for possibilities of going beyond the setup with fixed reference metric.

* The reference metric $f_{\mu\nu}$ in the general massive gravity action (4.2), or $\eta_{\mu\nu}$ in the dRGT theory (4.7), is not dynamically determined. It is put into the theory by hand and there is no obvious fundamental principle determining its form. From a field theoretical point of view, this and the related fact that the theory (without additional fields) breaks diffeomorphism invariance are rather undesirable features and it is not entirely clear how to interpret the presence of the fixed reference metric. The notion of a pre-geometric structure clearly goes against the main spirit of GR.

* Possibly related to the above point is the occurrence of superluminal and in particular acausal propagation in massive gravity<fnref n="21" /> [102, 122–126]. This raises serious questions for the physical viability of the theory, but see [2] for a discussion of some counter-arguments.

* The equations of motion of massive gravity (with flat or general reference metric) cannot give rise to homogeneous and isotropic solutions that lead to a viable cosmology [127–133]. Lacking this feature, the model cannot serve as a serious alternative to GR.

* More generally, in the parameter space of massive gravity there is no good limit which brings the equations and their solutions close to those of GR, which are well-tested. In principle, this issue could be resolved by the Vainshtein mechanism, but this causes serious tension with the cosmological Higuchi bounds [132, 134]. The underlying reason is that the metric that couples to matter contains additional degrees of freedom with respect to GR. This generically changes the physics significantly, as becomes apparent through the vDVZ discontinuity already at the linearised level.

From our point of view, the above shortcomings strongly motivate the extension of the ghost-free massive gravity action to a fully dynamical bimetric theory, which we shall focus on in the remainder of this article.

<footnote n="21">²¹ The analyses of these references do not extend to bimetric theory but in principle it is possible that similar problems arise there as well.</footnote>
## 5. Ghost-free bimetric theory

A natural question that arises within the formulation of massive gravity with general reference metric $f_{\mu\nu}$ is whether this second metric could be dynamical on its own without spoiling the consistency of the theory. The fact that one can extend the theory by a kinetic term and an equation of motion for $f_{\mu\nu}$ without reintroducing the Boulware–Deser ghost is not immediately evident from the consistency proofs of [59, 61, 63]. Nevertheless, Hassan and Rosen were able to show that one can indeed augment the action by an Einstein–Hilbert term for the second metric and let it be determined dynamically [76]. It has subsequently been confirmed that this also seems to be the unique kinetic term which can be added to the nonlinear massive gravity action in order to give dynamics$^{22}$ for $f_{\mu\nu}$ [135–137]. In this section we shall present the resulting ghost-free bimetric theory, review its consistency proof and discuss some of its most important features.

### 5.1. Action and equations of motion

The ghost-free action for Hassan–Rosen bimetric theory is given through [76]

$$S_{\text{HR}} = m_g^2 \int d^4 x \sqrt{g} R(g) + m_f^2 \int d^4 x \sqrt{f} R(f) - 2m_f^4 \int d^4 x \sqrt{g} \sum_{n=0}^4 \beta_n e_n (\sqrt{g^{-1}f}), \quad (5.1)$$

where $m_f$ is the ‘Planck mass’ for $f_{\mu\nu}$. It corresponds to the massive gravity action (4.2), extended by an Einstein–Hilbert term for the reference metric $f_{\mu\nu}$ which is promoted to a dynamical field on the same footing as $g_{\mu\nu}$. Consistency of the action now of course requires the inverse $f^{\mu\nu}$ to exist which, in analogy with GR, we shall always assume except perhaps at isolated physical singularities. The corresponding equations of motion for the two metrics read

$$R_{\mu\nu}(g) - \frac{1}{2}g_{\mu\nu}R(g) + \frac{m_f^4}{m_g^2}V_{\mu\nu}^g(g, f; \beta_n) = 0, \quad (5.2a)$$

$$R_{\mu\nu}(f) - \frac{1}{2}f_{\mu\nu}R(f) + \frac{m_f^4}{m_g^2}V_{\mu\nu}^f(g, f; \beta_n) = 0, \quad (5.2b)$$

in which the contributions from the interaction potential in terms of $S = \sqrt{g^{-1}f}$ are of the form

$$V_{\mu\nu}^g(g, f; \beta_n) = g_{\mu\rho} \sum_{n=0}^{3} (-1)^n \beta_n (Y_{(n)})^\rho_\nu (S), \quad (5.3)$$

$$V_{\mu\nu}^f(g, f; \beta_n) = f_{\mu\rho} \sum_{n=0}^{3} (-1)^n \beta_{4-n} (Y_{(n)})^\rho_\nu (S^{-1}). \quad (5.4)$$

The matrix functions $Y_{(n)}(S)$ have already been defined in (4.5). Due to the overall covariance of the interaction potential there is an identity between the divergences of these interaction contributions [70] (note that this identity follows from covariance and is otherwise independent of the form of the bimetric interactions)

<footnote>$^{22}$ Note that uniqueness of course only holds up to field redefinitions. The precise statement should therefore be that, assuming an Einstein–Hilbert term for $g_{\mu\nu}$ and a potential of the form given in (4.2), the only possible kinetic term for $f_{\mu\nu}$ is also of the Einstein–Hilbert form. Here we do not discuss the possibility of including additional non spin-2 degrees of freedom. For example, by adding extra scalar degrees of freedom one can consider $f(R)$ extensions of the theory (see e.g. [138–142]).</footnote>
$$ \sqrt{g} g^{\mu\rho} \nabla_{\rho} V_{\mu\nu}^{g} = - \sqrt{f} f^{\mu\rho} \tilde{\nabla}_{\rho} V_{\mu\nu}^{f}, \qquad (5.5) $$

where $\tilde{\nabla}$ is the covariant derivative compatible with $f_{\mu\nu}$. Due to this identity there is only one set of independent Bianchi constraints just as in massive gravity. Two important remarks are in order: firstly, note that the $\beta_4$-term $\sqrt{g} e_4(S) = \sqrt{f}$ is no longer non-dynamical but now contributes to the $f_{\mu\nu}$ equations of motion. Secondly, due to the more general symmetry property of the elementary symmetric polynomials (which is just a rewriting of (A.8))

$$ \sqrt{g} e_n(\sqrt{g^{-1}f}) = \sqrt{f} e_{4-n}(\sqrt{f^{-1}g}), \qquad (5.6) $$

the structure of the above action is symmetric in the two metrics. At the level of equations, this symmetry manifests itself through the identity

$$ V_{\mu\nu}^{f}(g, f; \beta_n) = V_{\mu\nu}^{g}(f, g; \beta_{4-n}). \qquad (5.7) $$

The metrics are therefore treated on the same footing and in section 5.4 we will see how a 'physical metric' is selected only by the choice of matter couplings. Various further aspects of the Hassan-Rosen action will be discussed in more detail in the remaining parts of this review.

## 5.2. Absence of ghost

The interaction potential of bimetric theory breaks the two diffeomorphism invariances of $g_{\mu\nu}$ and $f_{\mu\nu}$ down to their diagonal subgroup. In other words, the bimetric action is not invariant under independent coordinate transformations of the two metrics but only under those $\delta x^\mu = \xi^\mu$ that transform both metrics simultaneously in the same way: $\delta_\xi g_{\mu\nu} = -2g_{\rho(\mu} \nabla_{\nu)} \xi^\rho$ and $\delta_\xi f_{\mu\nu} = -2f_{\rho(\mu} \tilde{\nabla}_{\nu)} \xi^\rho$, where $\tilde{\nabla}$ is again the covariant derivative compatible with $f_{\mu\nu}$. For a general interaction potential, the degree of freedom counting therefore goes as follows: there are $2 \times 10 = 20$ components to start with; $2 \times 4 = 8$ of these get removed by gauge constraints and gauge fixing. Just as in the massive gravity case there is one set of Bianchi constraints, which can be taken to be either $\nabla^\mu V_{\mu\nu}^g = 0$ or $\tilde{\nabla}^\mu V_{\mu\nu}^f = 0$, since these are related by the identity(5.5). These vector constraints thus removes four additional degrees of freedom, leaving us with a total number of eight propagating modes. These correspond to the two degrees of freedom of a massless spin-2 field, five of a massive spin-2 and one additional scalar which gives rise to the Boulware-Deser ghost instability. In a consistent bimetric theory we therefore also need an additional constraint that eliminates the ghost mode.

In the ADM language, this means that we need the action to be linear in both lapses $N$ and $L$ of the two metrics as well as one set of three-dimensional vectors. In total there will then be five non-dynamical variables whose equations of motion become constraints: four corresponding to the gauge constraints associated to the overall diffeomorphism invariance and one extra constraint that removes the Boulware-Deser ghost.

Using the same variables as for massive gravity (see section 3.4) resulting from the redefinition (3.19), Hassan and Rosen were able to show that the bimetric action (5.1) indeed assumes the form [76]

$$ S_{\text{HR}} = \int d^4 x (L_i C^i + L \tilde{C} + N C), \qquad (5.8) $$

in which the constraints $C^i$, $\tilde{C}$ and $C$ are independent of the lapses $N$ and $L$ and the shifts $L^i$. As in the massive gravity case, the action is nonlinear in the redefined shift vectors $n^i$, which is a consequence of the breaking of one set of GCTs. The scalar constraints $\tilde{C}$ and $C$ contain terms coming from the interaction potential, whereas the vector constraints $C^i$ entirely originate from the Einstein-Hilbert terms Note that due to the redefinition of the form
$N^i = Ln^i + L^i + ND^i_k n^k$, all constraints receive contributions from the Einstein-Hilbert term for $g_{\mu\nu}$, which originally contains a term $N_i R^i$, see (3.4).

These results have subsequently been confirmed by explicit calculations and independently verified in various other approaches (see e.g. [63, 67–69]) which we shall not review here. Instead we note that the key observation which motivated Hassan and Rosen to study the fully dynamical extension of massive gravity was that due to the symmetry property (5.6), linearity of the interaction potential in the lapse $N$ implies that it must also be linear<sup>23</sup> in $L$. The corresponding secondary constraint which removes the canonical momentum of the ghost mode was also shown to exist [62]. Bimetric theory therefore gives rise to the correct amount of constraints in order to propagate the $5 + 2 = 7$ degrees of freedom, corresponding to a massive and a massless spin-2 field.

At the nonlinear level there is no unique split of degrees of freedom into ‘massive’ and ‘massless’ ones. As discussed in the following, the mass eigenstates can only be properly defined in the linearised theory around particular backgrounds.

### 5.3. Mass spectrum

The notion of mass is intimately related to the isometries of space-time. Its definition is most concise in Minkowski space where mass arises as a Casimir invariant of the Poincaré isometry group and it is possible to generalise that concept to space-times with the same amount of symmetries, i.e. Anti-de Sitter and de Sitter isometries. For less symmetric spaces it becomes more difficult to obtain a clear identification of mass, but one option is to classify a field as massless or massive depending on its number of propagating degrees of freedom. More precisely, if there is a parameter which when taken to zero increases the amount of gauge redundancy and thus reduces the number of propagating degrees of freedom to that of a massless theory then that parameter can loosely be identified with mass. This notion is implicit when we use the term ‘nonlinear massive gravity’ or when we speak of ‘massless’ and ‘massive’ degrees of freedom in bimetric theory. Nevertheless, around their maximally symmetric background solutions, we expect the nonlinear theories to give rise to a well-defined mass spectrum. In bimetric theory, such backgrounds are most easily obtained by making a proportional ansatz for the metrics<sup>24</sup>.

#### 5.3.1. Proportional background solutions.

Probably the simplest and yet remarkably important class of solutions to the bimetric equations of motion in vacuum is obtained by making an ansatz that conformally relates the two metrics, $\bar{f}_{\mu\nu} = c(x)^2 \bar{g}_{\mu\nu}$, where $c(x)$ is a space-time dependent function. Having plugged this ansatz into the equations of motion, we first note that the Bianchi constraint, $\nabla^\mu V^g_{\mu\nu} = 0$, immediately enforces $c(x) = \text{const}$. This is simply because this equation reduces to a polynomial in $c$ with constant coefficients that multiplies $\partial_\nu c$. This restricts our ansatz to proportional metrics

$$ \bar{f}_{\mu\nu} = c^2 \bar{g}_{\mu\nu} \quad \text{with} \quad c = \text{const.} \qquad (5.9) $$

Using this in the bimetric equation (5.2) we find that they simply reduce to two copies of Einstein’s equations

<footnote n="23"><sup>23</sup> There could of course be terms proportional to $NL$, but a closer inspection of the structure of the matrices $A$ and $B$ in (3.27) as well as the interactions in (3.33) reveals that such terms do not arise.</footnote>

<footnote n="24"><sup>24</sup> For general parameter values, there can exist additional maximally symmetric solutions, see section 6.1, but since they do not correspond to proportional metrics, we do not expect the perturbations to be diagonalisable into spin-2 mass eigenstates.</footnote>
$$G_{\mu\nu}(\bar{g}) + \Lambda_g \bar{g}_{\mu\nu} = 0, \quad G_{\mu\nu}(\bar{g}) + \Lambda_f \bar{g}_{\mu\nu} = 0. \qquad (5.10)$$

In this we have defined the cosmological constant contributions arising from the interaction potential

$$\Lambda_g(c) \equiv \frac{m^4}{m_g^2}(\beta_0 + 3c\beta_1 + 3c^2\beta_2 + c^3\beta_3), \\ \Lambda_f(c) \equiv \frac{m^4}{(m_f c)^2}(c\beta_1 + 3c^2\beta_2 + 3c^3\beta_3 + c^4\beta_4) \qquad (5.11)$$

and the Einstein tensors are the same in both equations since the Einstein tensor is scale invariant, $G_{\mu\nu}(c^2\bar{g}) = G_{\mu\nu}(\bar{g})$. Proportional backgrounds thus simply correspond to solutions to Einstein's equations in GR. Importantly, it means that this class of solutions actually captures all the solutions of GR. The difference of the background equations (5.10) implies that

$$\Lambda_g(c) = \Lambda_f(c). \qquad (5.12)$$

From (5.11) it is clear that this equation is a polynomial in $c$ with coefficients depending on the $\beta_n$ parameters and $\alpha = m_f/m_g$. In general, it serves to determine the proportionality constant $c$ of our ansatz in terms of the parameters of the theory and thereby specifies the solution completely. An important exception to this generic situation is the partially massless case, which we discuss in section 7.

Apart from being able to capture all solutions of GR, the proportional background solutions are of particular interest because they allow for a definite mass spectrum of fluctuations around them. In general, the linear perturbation equations have a rather complicated structure because in order to derive them one needs to vary the square-root matrix $S = \sqrt{g^{-1}f}$. As we already discussed in section 4.3 for the massive gravity setup, this is always possible but, for backgrounds giving rise to matrices $\bar{g}^{-1}\bar{f}$ that do not commute with the fluctuations, the resulting expressions are lengthy. In particular, the equations will not contain a mass term with Fierz-Pauli structure which makes it difficult, if not impossible, to uniquely identify the massive field. In contrast, for the proportional backgrounds we have $\bar{g}^{-1}\bar{f} = c^2\mathbf{1}$ which does commute with any other matrix and hence drastically simplifies the perturbation equations which now will exhibit the Fierz-Pauli structure.

### 5.3.2. Spectrum of linear mass eigenstates.
We shall now consider small perturbations around the proportional backgrounds for both of the metrics

$$g_{\mu\nu} = \bar{g}_{\mu\nu} + \delta g_{\mu\nu}, \quad f_{\mu\nu} = c^2 \bar{g}_{\mu\nu} + \delta f_{\mu\nu}. \qquad (5.13)$$

The variation of the square-root matrix is easily obtained as

$$\bar{g}_{\mu\rho} \delta S_{\nu}^{\rho} = \frac{1}{2c} (\delta f_{\mu\nu} - c^2 \delta g_{\mu\nu}). \qquad (5.14)$$

Plugging these into the bimetric equations (5.2), keeping only terms linear in the fluctuations and using the background equations, we obtain

$$\bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma} \delta g_{\rho\sigma} - \bar{\Lambda}_g (\delta g_{\mu\nu} - \frac{1}{2} \bar{G}_{\mu\nu} \bar{G}^{\rho\sigma} \delta g_{\rho\sigma}) - N \bar{G}_{\mu\rho} (\delta S^{\rho}_{\nu} - \delta^{\rho}_{\nu} \delta S^{\alpha}_{\alpha}) = 0, \qquad (5.15a)$$

$$\bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma} \delta f_{\rho\sigma} - \bar{\Lambda}_g (\delta f_{\mu\nu} - \frac{1}{2} \bar{G}_{\mu\nu} \bar{G}^{\rho\sigma} \delta f_{\rho\sigma}) + \alpha^{-2} N \bar{G}_{\mu\rho} (\delta S^{\rho}_{\nu} - \delta^{\rho}_{\nu} \delta S^{\alpha}_{\alpha}) = 0. \qquad (5.15b)$$

Here $N$ depends on $\alpha$, $c$ and the $\beta_n$ parameters and the explicit dependence can be read off from the Fierz-Pauli mass below. Here, we have made use of the fact that the background
metric can be conveniently rescaled by a constant without changing the structure of the linearised equations and expressed the equations with respect to a redefined background metric

$$ \bar{G}_{\mu\nu} \equiv (1 + \alpha^2 c^2) \bar{g}_{\mu\nu}. \qquad (516) $$

We have also redefined the cosmological constant with respect to this background

$$ \bar{\Lambda}_g \equiv (1 + \alpha^2 c^2)^{-1} \Lambda_g, \qquad R_{\mu\nu}(\bar{G}) = \bar{\Lambda}_g \bar{G}_{\mu\nu} \qquad (517) $$

and expressed the kinetic operator $\bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma}$ with respect to the $\bar{G}_{\mu\nu}$ background

$$ \bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma} \delta G_{\rho\sigma} = -\frac{1}{2} [\delta_{\mu}^{\rho} \delta_{\nu}^{\sigma} \bar{\nabla}^2 + \bar{G}^{\rho\sigma} \bar{\nabla}_{\mu} \bar{\nabla}_{\nu} - \delta_{\mu}^{\rho} \bar{\nabla}^{\sigma} \bar{\nabla}_{\nu} - \delta_{\nu}^{\rho} \bar{\nabla}^{\sigma} \bar{\nabla}_{\mu} - \bar{G}_{\mu\nu} \bar{G}^{\rho\sigma} \bar{\nabla}^2 + \bar{G}_{\mu\nu} \bar{\nabla}^{\rho} \bar{\nabla}^{\sigma}] \delta G_{\rho\sigma}, \qquad (518) $$

in which $\bar{\nabla}$ is the covariant derivative compatible with $\bar{G}_{\mu\nu}$.

The above equations can now easily be diagonalised into an equation for a massless and a massive perturbation. To this end, consider the following combinations of metric fluctuations

$$ \delta G_{\mu\nu} \equiv \delta g_{\mu\nu} + \alpha^2 \delta f_{\mu\nu}, \qquad \delta M_{\mu\nu} \equiv \frac{1}{2c} (\delta f_{\mu\nu} - c^2 \delta g_{\mu\nu}). \qquad (5.19) $$

Now, taking the appropriate linear combinations of the original fluctuation equations for $\delta g_{\mu\nu}$ and $\delta f_{\mu\nu}$ decouples the massless from the massive mode. The resulting equations read [115, 143]

$$ \bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma} \delta G_{\rho\sigma} - \bar{\Lambda}_g (\delta G_{\mu\nu} - \frac{1}{2} \bar{G}_{\mu\nu} \delta G) = 0, \qquad (5.20a) $$

$$ \bar{\mathcal{E}}_{\mu\nu}^{\rho\sigma} \delta M_{\rho\sigma} - \bar{\Lambda}_g (\delta M_{\mu\nu} - \frac{1}{2} \bar{G}_{\mu\nu} \delta M) + \frac{\bar{m}_{\text{FP}}^2}{2} (\delta M_{\mu\nu} - \bar{g}_{\mu\nu} \delta M) = 0, \qquad (5.20b) $$

where $\delta G = \bar{G}^{\mu\nu}\delta G_{\mu\nu}$ and $\delta M = \bar{G}^{\mu\nu}\delta M_{\mu\nu}$. The Fierz-Pauli mass in these equations is given by

$$ \bar{m}_{\text{FP}}^2 = \frac{(1 + (\alpha c)^2)^2}{\alpha^2 c} N = \frac{m^4}{m_g^2} \frac{1}{\alpha^2 c^2} (c\beta_1 + 2c^2\beta_2 + c^3\beta_3). \qquad (5.21) $$

We remind that in all of the above expressions c is to be regarded as a function of the Planck masses and the $\beta_n$ parameters, determined by the background equation $\Lambda_g = \Lambda_f$.

As advertised, (5.20a) and (5.20b), respectively, describe a massless and a massive spin-2 fluctuation in (Anti-) de Sitter background. At the linear level around proportional backgrounds, one can therefore assign two degrees of freedom to a massless fluctuation $\delta G_{\mu\nu}$ and the remaining five to a massive fluctuation $\delta M_{\mu\nu}$. The linearised action in terms of the mass eigenstates is [143]

$$ S_{\text{PB}} = \frac{1}{2} \int d^4 x \left[ \delta G_{\mu\nu} \bar{\mathcal{E}}^{\mu\nu\rho\sigma} \delta G_{\rho\sigma} - \bar{\Lambda}_g \left( \delta G^{\mu\nu} \delta G_{\mu\nu} - \frac{1}{2} \delta G^2 \right) + \delta M_{\mu\nu} \bar{\mathcal{E}}^{\mu\nu\rho\sigma} \delta M_{\rho\sigma} \right. \\ \left. - \bar{\Lambda}_g \left( \delta M^{\mu\nu} \delta M_{\mu\nu} - \frac{1}{2} \delta M^2 \right) + \frac{\bar{m}_{\text{FP}}^2}{2} \left( \delta M^{\mu\nu} \delta M_{\mu\nu} - \bar{g}_{\mu\nu} \delta M^2 \right) \right]. \qquad (5.22) $$

The main reason for choosing the new background $\bar{G}_{\mu\nu}$ was to render the final action and equations as simple as possible. Alternatively, we could have written all of the above expressions with respect to the original background $\bar{g}_{\mu\nu}$, in terms of $\Lambda_g$ in (5.11) and a properly rescaled Fierz-Pauli mass, $m_{\text{FP}}^2 \equiv (1 + \alpha^2 c^2) \bar{m}_{\text{FP}}^2$.
## 5.4. Couplings to matter

So far we have been dealing with theories involving only spin-2 degrees of freedom. In order to be accessible to any type of observations or experiments, the fields need to interact with ordinary matter. It may not come as a surprise that the set of allowed matter couplings which do not reintroduce the Boulware–Deser ghost is very restricted.

### 5.4.1. Ghost-free matter couplings.
The only known couplings which can be added to the bimetric action (5.1) without exciting the ghost are

$$S_m = \int d^4x \sqrt{g} \mathcal{L}_m(g, \Phi_g) + \int d^4x \sqrt{f} \tilde{\mathcal{L}}_m(f, \Phi_f), \qquad (5.23)$$

where $\mathcal{L}_m$ and $\tilde{\mathcal{L}}_m$ are standard minimally coupled matter Lagrangians of the same form as in GR. $\Phi_g$ and $\Phi_f$ schematically stand for sets of matter fields of any kind. Importantly, it was shown in [144, 145] that it is not possible to couple the same (dynamical) matter field to both metrics using minimal couplings, and hence $\Phi_g$ and $\Phi_f$ must be entirely independent. That coupling a field to both $g_{\mu\nu}$ and $f_{\mu\nu}$ reintroduces the Boulware–Deser ghost can be understood as follows: in GR with only one metric, the matter action becomes linear in the lapse $N$ when written in the canonical variables for the metric and the matter fields. A simple calculation in an example with a free scalar field $\varphi$ shows that this happens because the variation of the action with respect to $\varphi$ and hence also the canonical momentum of the scalar depend on $N$. This $N$-dependence is such that the action in terms of canonical variables is linear in $N$. In the bimetric case, however, when the free scalar is coupled to both $g_{\mu\nu}$ and $f_{\mu\nu}$, its canonical momentum will depend on $N$ and $L$ in a more complicated way and the action will not become linear in both of the lapses. As a consequence, the constraint which removes the Boulware–Deser ghost is lost. Note that this argument does not exclude coupling pure interaction terms (without appearance of time derivatives) to both of the metrics. This possibility seems quite contrived, however, and we do not discuss it any further but instead focus on couplings to two entirely independent matter sectors.

In the equations of motion the matter couplings enter in the form of stress-energy tensors

$$T_{\mu\nu}^g = -\frac{1}{\sqrt{g}} \frac{\delta(\sqrt{g} \mathcal{L}_m(g, \Phi_g))}{\delta g^{\mu\nu}}, \qquad T_{\mu\nu}^f = -\frac{1}{\sqrt{f}} \frac{\delta(\sqrt{f} \mathcal{L}_m(f, \Phi_f))}{\delta f^{\mu\nu}}. \qquad (5.24)$$

The bimetric equations (5.2) in the presence of matter sources thus become

$$\mathcal{G}_{\mu\nu}(g) + \frac{m^4}{m_g^2} V_{\mu\nu}^g(g, f; \beta_n) = \frac{1}{m_g^2} T_{\mu\nu}^g, \qquad (5.25a)$$

$$\mathcal{G}_{\mu\nu}(f) + \frac{m^4}{m_f^2} V_{\mu\nu}^f(g, f; \beta_n) = \frac{1}{m_f^2} T_{\mu\nu}^f. \qquad (5.25b)$$

We will comment on other matter couplings that have been studied in the literature below but first we shall discuss two important limits in the parameter space of bimetric theory which take it close to either GR or nonlinear massive gravity.

### 5.4.2. GR limit.
As we pointed out in section 4.5, the parameter space of massive gravity with a fixed reference metric does not include any region which is obviously close to GR. In the following we will see that, in contrast, a well-defined GR limit does exist for the bimetric theory. Since the structure of the bimetric action is completely symmetric in $g_{\mu\nu}$ and $f_{\mu\nu}$, either of the two metrics can play the role of the physical metric whose solutions will become
similar to those of GR. We choose this metric to be $g_{\mu\nu}$ here but note that everything can analogously be derived with the roles of the metrics interchanged.

It is straightforward to verify that the contributions from the interaction potential $V$ to the bimetric equations of motion (5.25) satisfy the following$^{25}$ identity [146, 147]

$$ \sqrt{g} g^{\mu\rho} V_{\rho\nu}^{g} + \sqrt{f} f^{\mu\rho} V_{\rho\nu}^{f} - \sqrt{g} V \delta_{\nu}^{\mu} = 0. \qquad (5.26) $$

Making use of this observation, the equations can be combined to give

$$ g^{\mu\rho} \mathcal{G}_{\rho\nu}(g) + \alpha^2 \det(\sqrt{g^{-1}f}) f^{\mu\rho} \mathcal{G}_{\rho\nu}(f) + \frac{m^4}{m_g^2} V \delta_{\nu}^{\mu} = \frac{1}{m_g^2} (g^{\mu\rho} T_{\rho\nu}^g + f^{\mu\rho} T_{\rho\nu}^f). \qquad (5.27) $$

This particular set of equations has interesting implications on classical solutions in bimetric theory which we shall come back to in section 6. For our purposes here it suffices to consider the dependence on the parameter $\alpha \equiv m_f/m_g$ in (5.27). In particular we note that in the limit $\alpha \to 0$, the equations reduce to [146, 147]

$$ g^{\mu\rho} \mathcal{G}_{\rho\nu}(g) + \frac{m^4}{m_g^2} V \delta_{\nu}^{\mu} = \frac{1}{m_g^2} (g^{\mu\rho} T_{\rho\nu}^g + f^{\mu\rho} T_{\rho\nu}^f). \qquad (5.28) $$

When $T_{\mu\nu}^f = 0$, i.e. when there is no matter sector for the metric $f_{\mu\nu}$, we can use the covariant derivative compatible with $g_{\mu\nu}$ to take the divergence of the equations which for a covariantly conserved source then gives $V$ = constant on-shell. In this case, the equations thus reduce to Einstein's equations for the physical metric $g_{\mu\nu}$ with Planck mass $m_g$ and cosmological constant $\frac{m^4}{m_g^2}V$

$$ \mathcal{G}_{\mu\nu}(g) + \frac{m^4}{m_g^2} V g_{\mu\nu} = \frac{1}{m_g^2} T_{\mu\nu}^g. \qquad (5.29) $$

Hence, the GR limit of bimetric theory is defined by

$$ \alpha \equiv \frac{m_f}{m_g} \to 0, \quad m_g = \text{const.}, \quad T_{\mu\nu}^f = 0, \qquad (5.30) $$

in which case the solutions for the physical metric $g_{\mu\nu}$ coincide with those of GR$^{26}$. Interestingly, a large value for the physical Planck mass $m_g$ automatically implies that bimetric theory is close to its GR limit, provided that $m_f$ is of reasonable size compared to other relevant scales (such as the electroweak scale, for instance).

The effect of taking the above limit on the $f_{\mu\nu}$ equation (5.25b) is that it becomes purely algebraical, $V_{\mu\nu}^f = 0$. The generic solutions to this equation are proportional backgrounds $f_{\mu\nu} = c^2 g_{\mu\nu}$ with $c$ determined by the condition $\Lambda_f(c) = 0$, where the function $\Lambda_f$ is the same as in (5.11). Then the cosmological constant in (5.29) is given by

$$ \frac{m^4}{m_g^2} V = \frac{m^4}{m_g^2} (\beta_0 + 4\beta_1 c + 6\beta_2 c^2 + 4\beta_3 c^3 + \beta_4 c^4). \qquad (5.31) $$

We further observe that, in the GR limit, the fluctuation of the physical metric $g_{\mu\nu}$ becomes massless as expected. This can be seen from (5.19), which for $\alpha \to 0$ gives $\delta G_{\mu\nu} \to \delta g_{\mu\nu}$. It is an important feature because, at least to our present knowledge, only the couplings of the

<footnote>$^{25}$ In fact, this identity can be derived for any form of covariant potential $V$ and does therefore not rely on the specific structure of the ghost-free bimetric action.</footnote>

<footnote>$^{26}$ One may be worried about occurrences of strong coupling for small values of $\alpha$. However, as has been discussed in detail in [148], the strong coupling scale of the massive spin-2 mode in fact grows with decreasing $\alpha$.</footnote>
original metrics to matter are ghost-free, whereas on the other hand the gravitational metric needs to behave like a massless spin-2 field for phenomenological reasons. In particular, since the massive mode decouples from the source, bimetric theory in the GR limit does not suffer from the vDVZ discontinuity and hence does not need to rely on the Vainshtein mechanism, a requirement which usually challenges the phenomenological viability of massive gravity models.

To summarise, in the limit of small $\alpha$, bimetric theory can be viewed as a smooth deformation of GR because the massive mode decouples in the limit. In particular, the gravitational metric satisfies Einstein's equations modified by a small correction and its fluctuations are dominated by the massless spin-2 mode. The dynamics of the massive field essentially decouple from the observable matter sector and its presence manifests itself only through the constant term in (5.29) which for a sufficiently small spin-2 mass could give rise to cosmic self-acceleration [148]. On the other hand, one can tune the $\beta_n$ parameters such that large contributions to the effective cosmological constant cancel (at the cost of giving up technical naturalness à la 't Hooft) and thus, even for a very large spin-2 mass, its observable effects could remain small.

Note also that, in this setup, the massive spin-2 field could potentially be a suitable dark matter candidate: it interacts with the Standard Model fields only very weakly, but couples non-minimally to the gravitational metric. A remarkable feature of this scenario would be that the closeness of the theory to GR goes hand in hand with the decoupling of the dark matter field, and both get related to the largeness of the Planck scale $m_g$. For related approaches in the context of bimetric theory, where not the spin-2 field itself is considered as the dark matter candidate but additional fields are invoked, see [149–152].

**5.4.3. Massive gravity limit.** Let us now consider the limit opposite to the one above, namely $\alpha \to \infty$. In this limit, bimetric theory reduces to massive gravity with a GR reference metric, i.e. $f_{\mu\nu}$ solves the standard sourced Einstein equations [153–156]. To see this, consider once more the bimetric equations of motion (5.25), this time along with the following scalings

$$ \alpha \equiv \frac{m_f}{m_g} \rightarrow \infty, \qquad m_g = \text{const.}, \qquad \frac{1}{M^2} \tilde{T}_{\mu\nu} \equiv \frac{1}{m_f^2} T_{\mu\nu}^f = \text{const.}, \\ \beta_4' \equiv \frac{m_g^2}{m_f^2} \beta_4 = \text{const.}, \qquad \beta_n = \text{const.} \quad \text{for } n \le 3. \qquad (5.32) $$

Note that we introduced a new mass scale $M$ and stress-energy $\tilde{T}$ here because now we also require a scaling of the matter fields in the $f_{\mu\nu}$ sector in order to be able to treat solutions for $f_{\mu\nu}$ that solve Einstein's equations in the presence of matter (for more on this see e.g. discussions in [115, 146, 154]). Similarly, the scaling of $\beta_4$ is required to keep a cosmological constant term for $f_{\mu\nu}$. The remaining $\beta_n$ are not allowed to scale since this would destroy interactions in the $g_{\mu\nu}$ equations (recall that $\beta_4$ only appear in the $f_{\mu\nu}$ equations). The $g_{\mu\nu}$ equations (5.25a) are then unaffected but in the $f_{\mu\nu}$ equation (5.25b) all bimetric interaction terms drop out. Thus, $f_{\mu\nu}$ is now determined by

$$ G_{\mu\nu}(f) + \frac{m^4}{m_g^2} \beta_4' f_{\mu\nu} = \frac{1}{M^2} \tilde{T}_{\mu\nu}, \qquad (5.33) $$

which is an Einstein equation with cosmological constant $\frac{m^4}{m_g^2} \beta_4'$ and Planck mass $M$. The exact limit thus decouples the dynamics of $f_{\mu\nu}$ which is now determined by equations that do not involve $g_{\mu\nu}$. We can obtain a solution to these equations and use it to replace $f_{\mu\nu}$ in the $g_{\mu\nu}$
equations. Effectively, the $g_{\mu\nu}$ equations are therefore the same as the ones obtained from varying the massive gravity action (4.2), in which $f_{\mu\nu}$ is taken to be a fixed reference metric that solves Einstein's equations (5.33). This picture explains the emergence of the fixed reference metric in massive gravity since the latter can be viewed as a particular point in the parameter space of the fully dynamical bimetric theory.

From the existence of the above limit we can infer that the solution space of bimetric theory is richer than that of massive gravity. Solutions $(g, f)$ to the bimetric equations of motion can be of the form $(g' + \mathcal{O}(\alpha^{-1}), f' + \mathcal{O}(\alpha^{-1}))$, with $(g', f')$ taken to be $\alpha$-independent. Only for such solutions is the limit $\alpha \to \infty$ well-defined and result in the massive gravity configurations $(g', f')$. Other bimetric solutions, however, do not possess a well-defined massive gravity limit, i.e. they become singular for $\alpha \to \infty$. Such configurations for the metrics are known to exist and have no massive gravity counterpart, but rather constitute a distinct feature of bimetric theory. Simple examples of proportional solutions that become singular in the limit have been found in [147]. Moreover, see [157] for a study of Hawking-Moss instanton solutions in bimetric theory which do not seem to allow for a well-defined massive gravity limit.

To summarise, all solutions of massive gravity can be viewed as arisen from bimetric theory, whereas it is quite easy to find solutions of the bimetric theory which does not have any massive gravity counterpart. Therefore, since the limit is singular, care has to be taken when arguing that results which hold in the massive gravity limit also hold in the full theory.

It has occasionally been argued that the above limiting procedure is somehow inferior and that the limit should instead be taken in the action. But, in fact, the two procedures of taking the limit in the equations or in the action are fully compatible when treated correctly. To see this let us briefly review the limiting procedure at the level of the action, as discussed, for instance, in [2]. In this case one starts by expanding the bimetric action (5.1) around, e.g., a constant curvature solution<fnref n="27" />

$$f_{\mu\nu} \rightarrow f_{\mu\nu}^{E} + \frac{\delta f_{\mu\nu}}{m_f} \qquad (5.34)$$

and then considers the limit $m_f \to \infty$. Since $f_{\mu\nu}^E$ is assumed to be a constant curvature metric the kinetic term for $f_{\mu\nu}$ simply reduces to the quadratic action for canonically normalised fluctuations on this background, i.e.

$$m_f^2 \sqrt{f} R(f) \rightarrow \frac{1}{2} \delta f_{\mu\nu} \tilde{\mathcal{E}}^{\mu\nu\rho\sigma} \delta f_{\rho\sigma}, \qquad (5.35)$$

where $\tilde{\mathcal{E}}^{\mu\nu\rho\sigma}$ is the operator defined in (2.12), written in terms of the metric $f_{\mu\nu}^E$ and its curvatures. In the interaction potential on the other hand, all the $\delta f_{\mu\nu}$ fluctuations decouple in the limit and, assuming an appropriate scaling of $\beta_4$, only a cosmological term for the fluctuations remains. One then ends up with a non-covariant action in terms of a decoupled linear spin-2 field $\delta f_{\mu\nu}$ and a nonlinear spin-2 field $g_{\mu\nu}$ whose interaction potential contains the fixed reference metric $f_{\mu\nu}^E$. Now varying with respect to the dynamical fields $g_{\mu\nu}$ and $\delta f_{\mu\nu}$ results in the massive gravity equations for $g_{\mu\nu}$ (containing $g_{\mu\nu}$ and $f_{\mu\nu}^E$), supplemented with a completely decoupled linear equation for $\delta f_{\mu\nu}$. Of course, this procedure is only self-consistent provided that $f_{\mu\nu}^E$ is really a finite constant curvature background solution in the limit. Ensuring this leads exactly to the limiting procedure of the equations we have discussed above. It is also straightforward to see that if we express our equations (5.25a) and (5.33) in

<footnote n="27">²⁷ More generally one can couple appropriately scaled matter fields to the $f_{\mu\nu}$ metric and consider expanding around a solution of (5.33), but this does not change the main arguments.</footnote>
the limit in terms of $g_{\mu\nu}$, $\delta f_{\mu\nu}$ and $f_{\mu\nu}^E$, then we end up with exactly the same result as that obtained from the action (the reason being that, if $\delta f_{\mu\nu}/(\alpha m_g)$ is assumed subdominant to $\delta g_{\mu\nu}/m_g$ in the action, the same will of course hold in the equations). Therefore taking the limit in the equations is equivalent to doing it in the action, but the former procedure deals more directly with solutions and the requirements for these to exist.

**5.4.4. Other matter couplings.** Matter couplings differing from the consistent ones in (5.23) have also been studied in the literature. Particular attention has been paid to the 'doubly coupled' theory for which the matter sectors of $g_{\mu\nu}$ and $f_{\mu\nu}$ contain the same fields [158, 159] but, as explained above, these couplings reintroduce the Boulware–Deser ghost [144, 145]. On the other hand, it turns out that a certain combination of the two metrics can be coupled to matter without exciting the ghost in the low-energy effective field theory [145]. This 'effective metric' contains two arbitrary parameters a and b and is of the form

$$G_{\mu\nu}^{\text{eff}} = a^2 g_{\mu\nu} + 2ab g_{\mu\rho} (\sqrt{g^{-1}f})^{\rho}_{\nu} + b^2 f_{\mu\nu}, \qquad (5.36)$$

whose 'uniqueness' has been discussed<fnref n="28" /> in [161, 162]. The reason this metric is special is because it can be written, in matrix notation, $G^{\text{eff}} = a^2 g(1 + (b/a)S)^2$, and hence $\sqrt{\det(G^{\text{eff}})} = a^4 \sqrt{g} \det(1 + (b/a)S)$. The identity (A.5) then reveals that any vacuum contribution generated by matter coupled to such a $G^{\text{eff}}$ will not alter the ghost-free form of the bimetric interaction potential. Interestingly, this metric can be written as a Finsler metric [159] and therefore provides a situation where the geometric interpretation is shifted from the standard (pseudo) Riemannian description to an effective Finsler geometry<fnref n="29" />. Since our main interest here lies in working with the full bimetric action and not in the effective field theory picture, we will not discuss phenomenological implications of such a matter coupling in this review. The interested reader is referred to the large variety of [99, 151, 152, 165–174] in which the effective coupling has been studied further, mostly in the context of cosmology.

## 6. Classical solutions

Despite the fact that the consistent theories have only been known for a few years, there already exists an extensive literature on classical solutions in ghost-free massive gravity and bimetric theory. Here we focus on bimetric theory whose solution spectrum is richer than that of massive gravity (see section 5.4.3). Perhaps most interesting is the class of spherically symmetric solutions, which can potentially be used to study stars, galaxies, black holes and cosmology. An immediate problem that arises in this context is that the bimetric theory has no known analogue of Birkhoff's theorem<fnref n="30" /> and therefore many of the valuable uniqueness theorems of GR do not apply. This means that in many situations several solutions may exist for the same problem and one is left with the task of sorting out the most relevant one. A strong physically motivated guide here is to explore the stability of solutions under perturbations. Another complication is the analytical complexity of the nonlinear equations of motion, and in many cases only numerical solutions are known as of yet. In this review we have mainly kept our attention towards analytically tractable problems and will continue this

<footnote n="28">²⁸ More precisely, the metric (5.36) is unique up to multiplication of the right-hand side by an arbitrary matrix of unit determinant. When it is expressed in terms of vierbeins more ambiguities arise [160].</footnote>

<footnote n="29">²⁹ A Finsler geometry departs from the pseudo Riemannian geometry in that it characterises a manifold with a norm but it is not necessarily infinitesimally Minkowski [163, 164].</footnote>

<footnote n="30">³⁰ Which implies absence of monopole radiation in GR. See e.g. [175] for a nice GR oriented discussion which also gives appropriate credit to earlier independent findings of this important theorem.</footnote>
practice. We will therefore restrict our discussion mainly to black hole and cosmological solutions, discuss features of these which can clearly be discerned analytically and only comment briefly on some phenomenological issues. For the spherically symmetric solutions and in particular their applications to black hole studies there already exist a few good reviews on the current status [176–178]. We refer the interested reader to these for additional details. Before discussing particular features of the spherically symmetric solutions, we make some general remarks.

## 6.1. General properties

Recall the form of the bimetric equations of motion in the presence of matter sources

$$ \mathcal{G}_{\mu\nu}(g) + \frac{m^4}{m_g^2} V_{\mu\nu}^g(g, f; \beta_n) = \frac{1}{m_g^2} T_{\mu\nu}^g, \qquad (6.1a) $$

$$ \mathcal{G}_{\mu\nu}(f) + \frac{m^4}{m_f^2} V_{\mu\nu}^f(g, f; \beta_n) = \frac{1}{m_f^2} T_{\mu\nu}^f, \qquad (6.1b) $$

with the contributions from the interaction potential given in (5.3). It is clear that due to the presence of these additional interaction terms, in general, the solutions to the bimetric equations will significantly differ from those obtained in GR. From a phenomenological point of view, such large new effects are not desirable since Einstein's theory is well-tested over a wide range of distance scales. In order not to be ruled out immediately, any modification of gravity needs to give rise to solutions that do not deviate too much from those of GR. As we already saw in section 5.4.2, in the limit of small $\alpha \equiv m_f/m_g$ and a vanishing source for $f_{\mu\nu}$, all equations for the physical metric $g_{\mu\nu}$ will smoothly approach the form of the GR equations and so will their linear perturbations. This is already good news for the viability of bimetric theory. On the other hand, it may also render the theory less interesting because if all solutions are (almost) indistinguishable from those of GR in the weak-field limit then there is little room for predicting new signatures that could be observed. It is therefore interesting to see if the above equations, away from the GR limit, can still give rise to solutions that behave similarly to those of Einstein's theory.

In section 5.3.1 we already encountered the proportional backgrounds as an example of GR solutions in bimetric theory. However, in the presence of matter, these backgrounds only exist if the stress-energy tensors satisfy the rather strict constraint $T_{\mu\nu}^f = \alpha^2 T_{\mu\nu}^g$. As we will see below, the bimetric equations can also reduce to Einstein equations after inserting ansätze for the metrics which possess particular symmetry properties (e.g. spherical symmetry). In these cases, although the background solutions reproduce exactly the predictions of GR, differences will generically occur at the level of perturbations and the theory can in principle make testable predictions.

In this context, it is important to notice one more feature of bimetric theory: if either $g_{\mu\nu}$ or $f_{\mu\nu}$ is assumed to be an exact solution to the Einstein equations, then the bimetric equations will force the second metric to also solve (a different set of) Einstein's equations<fnref n="31" /> [147]. In vacuum this can be seen as follows: recall that using the identity (5.26), we were able to combine the bimetric equations into

<footnote n="31">³¹ In fact, this statement does not depend on the ghost-free structure but holds for any covariant bimetric interaction potential.</footnote>
$$g^{\mu\rho}G_{\rho\nu}(g) + \alpha^2 \det(\sqrt{g^{-1}}f) f^{\mu\rho}G_{\rho\nu}(f) + \frac{m^4}{m_g^2}V\delta^\mu_\nu = 0. \quad (6.2)$$

If in this equation we assume that, for instance, $f_{\mu\nu}$ is a GR solution with cosmological constant $\tilde{\Lambda}$

$$G_{\mu\nu}(f) = -\tilde{\Lambda}f_{\mu\nu}, \quad (6.3)$$

then (6.2) takes the form

$$g^{\mu\rho}G_{\rho\nu}(g) - \alpha^2 \det(\sqrt{g^{-1}}f)\tilde{\Lambda}\delta^\mu_\nu + \frac{m^4}{m_g^2}V\delta^\mu_\nu = 0. \quad (6.4)$$

Taking the covariant divergence of this equation implies that the terms proportional to $\delta^\mu_\nu$ are constant and hence the equation reduces to an Einstein equation for $g_{\mu\nu}$

$$G_{\mu\nu}(g) = -\Lambda g_{\mu\nu}, \quad \Lambda \equiv -\alpha^2 \det(\sqrt{g^{-1}}f)\tilde{\Lambda} + \frac{m^4}{m_g^2}V. \quad (6.5)$$

This proof can straightforwardly be generalised to the equations including matter sources. However, in this case GR solutions do not exist unless the sources for $g_{\mu\nu}$ and $f_{\mu\nu}$ are related in a particular way [147]. Solutions of this type (with stress-energy tensors assumed to resemble perfect fluids) have been found in [179].

Another interesting general feature of the bimetric solutions was noticed in [180]. There it was found that the effective stress-energy tensors that the bimetric interactions generate, i.e. the negative of $V_{\mu\nu}^g$ and $V_{\mu\nu}^f$, never satisfy the null-energy condition simultaneously unless the two metrics are proportional. In more detail [180], found that if $k^\mu$ is a null-vector with respect to $g_{\mu\nu}$, i.e. $g_{\mu\nu}k^\mu k^\nu = 0$, then one can always define $\tilde{k}^\mu = (S^{-1})^\mu_\rho k^\rho$ which turns out to be a null-vector with respect to $f_{\mu\nu}$, i.e. $f_{\mu\nu}\tilde{k}^\mu\tilde{k}^\nu = g_{\mu\nu}k^\mu k^\nu = 0$. This follows since $f_{\mu\nu} = g_{\mu\rho}S^\rho_\sigma S^\sigma_\nu = g_{\sigma\rho}S^\rho_\mu S^\sigma_\nu$ and can be used to demonstrate that if e.g. $V_{\mu\nu}^g$ satisfy the null-energy condition $V_{\mu\nu}^g k^\mu k^\nu \le 0$ (where the direction of the inequality follows from our sign convention in defining $V_{\mu\nu}^g$) then $V_{\mu\nu}^f$ will satisfy an opposite inequality $V_{\mu\nu}^f \tilde{k}^\mu \tilde{k}^\nu \ge 0$. Furthermore the inequalities only saturate for proportional solutions, when the interactions reduce to pure cosmological constant contributions (see section 5.3.1). Of course the null-energy condition is usually applied to the matter sector within the context of GR so it is not completely obvious what a violation actually means physically in this case. For example, if we consider the case where only $g_{\mu\nu}$ couples to a matter source $T_{\mu\nu}^g$ and take it that $V_{\mu\nu}^f$ does not (does) satisfy the null-energy condition. Then generically $V_{\mu\nu}^g$ will (will not) satisfy it which implies that the standard physical interpretation of the null-energy condition on the source term $T_{\mu\nu}^g$ may change. This interesting observation is something which deserves further study. One immediate consequence however was the prediction of worm-hole solutions in the theory, which were subsequently found and analysed in [181].

## 6.2. Black hole solutions

Spherically symmetric solutions were first discussed in [182, 183] and have since then received considerable attention [184–193]. As already mentioned in the introduction of this section, there is no known analogue of Birkhoff’s theorem in bimetric theory and therefore many of the uniqueness theorems of GR fail straight away. In fact, for spherically symmetric ansätze, it is fairly straightforward to make an initial separation into two classes of general solutions, the so-called bidiagonal and non-bidiagonal solutions. These labels are not that
imaginative but they do keep to the point: the bidiagonal solutions are solutions for which both metrics can be brought to a simultaneously diagonal form while the non-bidiagonal solutions are those solutions which cannot. Another important class of solutions are the proportional solutions discussed in section 5.3.1. Now, we can always use the isometries to choose coordinates such that at least one of the metrics is diagonal and therefore the proportional solutions fall into the broader class of bidiagonal solutions.

The most general spherically symmetric ansätze can be written in the form

$$ 
\begin{aligned}
g_{\mu\nu}dx^\mu dx^\nu &= -U^2dt^2 + V^2dr^2 + r^2d\Omega^2, \\
f_{\mu\nu}dx^\mu dx^\nu &= -A^2dt^2 + B^2dr^2 + C^2dtdr + D^2d\Omega^2,
\end{aligned}
\qquad (6.6) $$

where $U, V, A, B, C, D$ are all functions of the radial ($r$) and temporal ($t$) coordinates. We have made use of the spherical isometry to fix the angular and radial coordinates such that the angular measure $d\Omega^2 = d\theta^2 + \sin^2\theta d\phi^2$ is in the standard form and comes with the standard factor of $r^2$ for $g_{\mu\nu}$, which also puts $g_{\mu\nu}$ in diagonal form. This form has been heavily used in the literature and serves to simplify parts of our discussions, but it should be noted that in some situations, concerning in particular black hole solutions, it is preferable to work with coordinates that are regular at the horizon, such as the advanced Eddington-Finkelstein coordinates (see e.g. [178, 190]). Since $g_{\mu\nu}$ is diagonal it follows from the equations of motion

$$ g^{\rho\mu}G_{\mu\nu} + \frac{m^4}{m_g^2}g^{\rho\mu}V_{\mu\nu}^g = 0, \qquad (6.7) $$

that $g^{\rho\mu}V_{\mu\nu}^g$ must be diagonal on the solution. The only off-diagonal terms in $g^{\rho\mu}V_{\mu\nu}^g$ that do not vanish identically are

$$ g^{\mu\nu}V_{\mu\nu}^g \sim g^{\mu\nu}V_{\mu t}^g \sim C \left( \beta_1 + 2\beta_2 \frac{D}{r} + \beta_3 \frac{D^2}{r^2} \right). \qquad (6.8) $$

This can vanish either if $C = 0$, which characterises the bidiagonal solutions, or if $D = \omega r$ with $\omega$ being a solution of

$$ \beta_1 + 2\beta_2\omega + \beta_3\omega^2 = 0. \qquad (6.9) $$

Notice that the second option, which characterises the non-bidiagonal solutions, does not exist if two of the $\beta_n$ parameters with $n = 1, 2, 3$ are zero. Furthermore, the condition forces $\omega$ to be a constant. We will now discuss the two classes of solutions in some more detail.

### 6.2.1. Bidiagonal black hole solutions.
We start by discussing the bidiagonal solutions, which are obtained from (6.6) with $C = 0$. As a simple but interesting and illustrative warmup we recall the proportional background solutions of section 5.3.1, with $f_{\mu\nu} = c^2g_{\mu\nu}$. If either metric is diagonal, then these are obviously bidiagonal. It should also be clear that any vacuum solutions of the standard GR equations will in fact also be solutions of the bimetric equations for the proportional backgrounds (with appropriate restrictions on the constant of proportionality). In particular we may now consider bidiagonal black hole metrics. Since these are proportional vacuum solutions the linearisation around these solutions follows the analysis of section 5.3.2 and the fluctuations obey the massless equations

$$ \delta E_{\mu\nu}^0 = \mathcal{E}_{\mu\nu}^{\rho\sigma} \delta G_{\rho\sigma} - \Lambda \left( \delta G_{\mu\nu} - \frac{1}{2} g_{\mu\nu} g^{\rho\sigma} \delta G_{\rho\sigma} \right) = 0 \qquad (6.10) $$
and the massive equations

$$ \delta E_{\mu\nu}^{M} = \mathcal{E}_{\mu\nu}^{\rho\sigma} \delta M_{\rho\sigma} - \Lambda \left( \delta M_{\mu\nu} - \frac{1}{2} g_{\mu\nu} g^{\rho\sigma} \delta M_{\rho\sigma} \right) + \frac{m_{\text{FP}}^2}{2} (\delta M_{\mu\nu} - g_{\mu\nu} g^{\rho\sigma} \delta M_{\rho\sigma}) = 0, \qquad (6.11) $$

where, in terms of the original fluctuations $\delta g_{\mu\nu}$ and $\delta f_{\mu\nu}$, the massless field is given by $\delta G_{\mu\nu} = \delta g_{\mu\nu} + \alpha^2 \delta f_{\mu\nu}$, and the massive field is obtained from $2c \delta M_{\mu\nu} = \delta f_{\mu\nu} - c^2 \delta g_{\mu\nu}$. We also recall that (with $\nabla_\mu$ associated to the background field $g_{\mu\nu}$)

$$ \mathcal{E}_{\mu\nu}^{\rho\sigma} \delta M_{\rho\sigma} = -\frac{1}{2} [\delta_{\mu}^{\rho} \delta_{\nu}^{\sigma} \nabla^2 + g^{\rho\sigma} \nabla_{\mu} \nabla_{\nu} - \delta_{\mu}^{\rho} \nabla^{\sigma} \nabla_{\nu} - \delta_{\nu}^{\rho} \nabla^{\sigma} \nabla_{\mu} \\ - g_{\mu\nu} g^{\rho\sigma} \nabla^2 + g_{\mu\nu} \nabla^{\mu} \nabla^{\sigma}] \delta M_{\rho\sigma}. \qquad (6.12) $$

The massless equations (6.10) are precisely those of linearised GR and a standard treatment shows that they propagate the two ($\pm 2$) helicity states of a massless spin-2 field. Let us therefore focus on the massive field equations (6.11). Due to the Bianchi identities obeyed by the Einstein tensor (including the $\Lambda$ terms) a covariant divergence of these yields

$$ \nabla^{\mu} \delta E_{\mu\nu}^{M} = \frac{m_{\text{FP}}^2}{2} (\nabla^{\mu} \delta M_{\mu\nu} - g^{\rho\sigma} \nabla_{\nu} \delta M_{\rho\sigma}), \qquad \nabla^{\mu} \delta M_{\mu\nu} = g^{\rho\sigma} \nabla_{\nu} \delta M_{\rho\sigma}, \qquad (6.13) $$

where the last equality follows on-shell and provides four non-dynamical constraint equations for $\delta M_{\mu\nu}$. A second covariant divergence of these equations yields

$$ \nabla^{\mu} \nabla^{\nu} \delta E_{\mu\nu}^{M} = \frac{m_{\text{FP}}^2}{2} (\nabla^{\mu} \nabla^{\nu} \delta M_{\mu\nu} - g^{\rho\sigma} \nabla^2 \delta M_{\rho\sigma}). \qquad (6.14) $$

Tracing instead the field equations (6.11) we get

$$ g^{\mu\nu} \delta E_{\mu\nu}^{M} = g^{\mu\nu} \nabla^2 \delta M_{\mu\nu} - \nabla^{\mu} \nabla^{\nu} \delta M_{\mu\nu} + \left( \Lambda - \frac{3m_{\text{FP}}^2}{2} \right) g^{\mu\nu} \delta M_{\mu\nu}. \qquad (6.15) $$

It then follows that the particular linear combination

$$ 2 \nabla^{\mu} \nabla^{\nu} \delta E_{\mu\nu}^{M} + m_{\text{FP}}^{2} g^{\mu\nu} \delta E_{\mu\nu}^{M} = \frac{m_{\text{FP}}^{2}}{2} (2\Lambda - 3m_{\text{FP}}^{2}) g^{\mu\nu} \delta M_{\mu\nu}, \qquad (6.16) $$

constitutes another on-shell scalar<fnref n="32" /> constraint, $g^{\mu\nu} \delta M_{\mu\nu} = 0$, which then also implies $\nabla^{\mu} \delta M_{\mu\nu} = 0$. Together these equations correspond to five constraints that can be used to remove five degrees of freedom from the original ten components of the symmetric tensor perturbations $\delta M_{\mu\nu}$. Enforcing these constraints, the equations of motion for the massive spin-2 fluctuations can therefore be reduced to the following system

$$ (\nabla^2 - m_{\text{FP}}^2) \delta M_{\mu\nu} + 2 R_{\mu}^{\rho} \delta M_{\rho\nu} = 0, \qquad \nabla^{\mu} \delta M_{\mu\nu} = 0, \qquad g^{\mu\nu} \delta M_{\mu\nu} = 0. \qquad (6.17) $$

This brief discussion provides a generalisation of the constraint analysis of section 2.2, when augmented to the case of constant curvature backgrounds, and also complements the analysis of section 5.3.2. Moreover, as first recognised in [195], the dynamical equations for $\delta M_{\mu\nu}$ in the above form are actually identical to the equations of a 5D black string when projected down onto its 4D sub-manifold, namely

<footnote n="32">³² An exceptional situation occurs when the parameters satisfy $2\Lambda = 3m_{\text{FP}}^2$, saturating the so called Higuchi bound [81]. In this case the left-hand side of (6.16) vanishes identically, such that a new linear scalar gauge symmetry emerges and the theory propagates a partially massless spin-2 field with only four degrees of freedom [85, 194]. We will get back to this in section 7.</footnote>
$$ (\nabla^2 - k^2)\delta M_{\mu\nu} + 2R_{\mu\nu}^{\rho\sigma} \delta M_{\rho\sigma} = 0, \qquad (6.18) $$

where $k^2$ denotes the wave number of the transverse Fourier mode. This equation is identical to (6.17) with the replacement $m_{FP} \to k$. It has immediate and interesting implications since it is well known that the solution of (6.18) exhibits the so called Gregory–Laflamme instability [196] (for a broad review on this subject see [197]). Namely, the solution is unstable in the parameter region $0 < k < k_c$, where $k_c$ is of the order of the inverse Schwarzschild radius. In the case of the 5D black string this instability concerns the mode longitudinal to the 5th direction along the black string and the end point of the instability is known to result in a sort of ‘pinching’ of the black string into an open necklace with beads of 4D Schwarzschild black holes strung on it. In the case of the bid diagonal bimetric solution when the metrics are proportional Schwarzschild (or Schwarzschild–de Sitter) metrics it follows that it too will have a similar instability$^{33}$. But since there is no 5th dimension, no such interpretation is available and the status of the end point is presently unknown (one possibility is mentioned below). It is clear however that the instability only manifests itself over a characteristic time-scale comparable to the Hubble time $H_0^{-1}$ for a small enough mass parameter $m_{FP} \sim H_0$ and it would therefore be hard to detect any signature connected with this. A similar study of the bid diagonal Kerr solutions have revealed that these additionally suffer from a super-radiant instability [187]. For more details and discussions on this we refer the interested reader to the reviews [176–178] and references therein.

Another simple and illustrative example is to consider the metric ansätze (6.6) as a static perturbation of proportional Minkowski$^{34}$ backgrounds, i.e. $\tilde{f}_{\mu\nu} = c^2\tilde{g}_{\mu\nu} = c^2\eta_{\mu\nu}$. Thus we consider the functions in (6.6), with $C = 0$ and no temporal dependence, as given by

$$ \begin{aligned} U &= 1 + \delta U, & V &= 1 + \delta V, \\ A &= c(1 + \delta A), & B &= c(1 + \delta B), & D &= c(r + \delta D). \end{aligned} \qquad (6.19) $$

The equations of motion for the perturbations can then be solved perturbatively to give (see e.g. [182, 184])

$$ \begin{aligned} \delta U &= \frac{M_1}{r} + \frac{M_2}{r} e^{-r m_{FP}}, \\ \delta V &= \frac{M_1}{r} - \frac{M_2(1 + r m_{FP})}{2r} e^{-r m_{FP}}, \\ \delta A &= \frac{M_1}{r} - \frac{M_2}{\alpha^2 c^2 r} e^{-r m_{FP}}, \\ \delta B &= \frac{M_1}{r} - (1 + r m_{FP})[\alpha^2 c^2 r^2 m_{FP}^2 + 2(1 + \alpha^2 c^2)] \frac{M_2}{2\alpha^2 c^2 m_{FP}^2 r^3} e^{-r m_{FP}}, \\ \delta D &= (1 + \alpha^2 c^2)[1 + r m_{FP} + r^2 m_{FP}^2] \frac{M_2}{2\alpha^2 c^2 m_{FP}^2 r^2} e^{-r m_{FP}}, \end{aligned} \qquad (6.20) $$

where $M_1$ and $M_2$ are integration constants and we recall that the Fierz–Pauli mass $m_{FP}$ is given by

<footnote>$^{33}$ An interesting exception is the Schwarzschild–de Sitter case with parameters fixed to the partially massless model discussed in section 7. For that case the instability is absent [185].</footnote>

<footnote>$^{34}$ Note that these backgrounds do not solve the bimetric equations in general, but the existence of such solutions requires fixing either $c$ or one of the $\beta_n$ parameters in terms of the others.</footnote>
$$m_{\text{FP}}^2 = \frac{m^4}{m_g^2} \frac{1 + \alpha^2 c^2}{\alpha^2 c^2} (c\beta_1 + 2c^2\beta_2 + c^3\beta_3). \qquad (6.21)$$

The above expressions in (6.20) clearly have the form of Yukawa suppressed GR solutions for spherically symmetric perturbations of flat space-time. They exhibit the vDVZ discontinuity since, in the small mass limit $m_{\text{FP}} \to 0$, the combination

$$\delta U + \delta V = \frac{M_2(1 - rm_{\text{FP}})}{2r} e^{-rm_{\text{FP}}} \to \frac{M_2}{2r}, \qquad (6.22)$$

does not vanish as it would in GR. On the other hand, we also see that both $\delta B$ and $\delta D$ diverge for $m_{\text{FP}} \to 0$, so this limit is in fact not well-defined on the solution. These divergences in the small-mass limit are remedied via the Vainshtein mechanism when going to higher orders in perturbation theory. Although the full nonlinear forms of these solutions are not known analytically their perturbative and numerical existence provides some initial hope of such a completion. The perturbative form, going up to second order, has been used for initial studies of e.g. strong lensing [184] and in confirming the Vainshtein mechanism in bimetric theory [13, 183, 192]. Making a full perturbative ansatz has also allowed to find numerical solutions to the nonlinear equations which are asymptotically flat and have massive hair. It has been conjectured that these solutions are the end point of the linear bidiagonal Schwarzschild instability discussed above [186]. Matching these solutions however requires the black hole to be of cosmological size and they are therefore unlikely to be of astrophysical interest.

In more general setups, away from proportional backgrounds, all known analytical solutions correspond to both metrics being of standard GR form and do not have massive hair. On the other hand, it is known that the bidiagonal class allows also for numerical solutions which are of non-GR form and do contain massive hair, typically with anti-de Sitter asymptotics. So far these more exotic solutions have only been studied numerically and we refer the reader to [177, 178] for more discussions of their explicit behaviour.

### 6.2.2. Non-bidiagonal black hole solutions.
The non-bidiagonal solutions have $C \neq 0$ and $D = \omega r$ in (6.6), with $\omega$ being a constant solution of (6.9). Furthermore, the Bianchi constraint $\nabla^\mu V_{\mu\nu}^g = 0$ implies that (see e.g. [177])

$$(\beta_2 + \beta_3\omega) \left[ (\omega - S_t^t)(\omega - S_r^r) + \left(\frac{U}{V} S_r^t\right)^2 \right] = 0, \qquad (6.23)$$

where $S_{\mu\nu}^\rho$ are components of the square-root matrix $S = \sqrt{g^{-1}f}$. This condition together with (6.9) suggests that the non-bidiagonal solutions are a very non-generic class in the sense that the metric coefficients are forced to have a quite peculiar dependence on the $\beta_n$ parameters of the action. If we assume that $\omega$ is a solution of (6.9) and that (6.23) is satisfied, we find that the equations of motion decouple into two separate Einstein equations

$$\mathcal{G}_{\mu\nu} + \frac{m^4}{m_g^2} \lambda_g g_{\mu\nu} = 0, \qquad \tilde{\mathcal{G}}_{\mu\nu} + \frac{m^2}{\alpha^2} \lambda_f f_{\mu\nu} = 0, \qquad (6.24)$$

where the cosmological terms are given by

$$\lambda_g = \beta_0 + 2\beta_1\omega + \beta_2\omega^2, \qquad \lambda_f = \frac{1}{\omega^2}(\beta_2 + 2\beta_3\omega + \beta_4\omega^2). \qquad (6.25)$$

Therefore, this peculiar way of fixing the metric functions in terms of the $\beta_n$ parameters is an alternative and more complicated way of tuning the interactions to be purely cosmological, as opposed to the more obvious proportional solutions. In this context, recall the results of
section 6.1, showing that whenever one metric is a solution to Einstein's equations, the other one must be an Einstein solution as well. That such solutions (both Einstein, but with different cosmological constants) exist could perhaps have been anticipated by the complicated matrix structure of the interaction terms All black hole solutions in the non-bidiagonal class are therefore of standard GR type. An initial treatment of perturbations around these solutions can be found in [191]. For radial modes the analysis simplifies drastically and it is possible to get analytical expressions for the perturbations. It turns out that, while these are regular at the horizon, they are not regular at infinity. This implies that unstable radial modes are not allowed on physical grounds and that the corresponding black holes may in fact be stable, at least in the perturbative sense. For a detailed discussion on the current status of perturbations of various known black hole solutions we refer the reader to [178].

### 6.3. Cosmological solutions

On small scales GR makes very good predictions and, from a phenomenological point of view, there is no need to look for a modification of the theory. On cosmological scales, however, it seems that either quantum field theory or gravity (or both) have to be modified in order to explain the observed value of the cosmic acceleration in a satisfactory way. In order for this review to be self-contained and to set up some notation we have provided a condensed summary of standard GR cosmology in appendix C.

Due to the conceptual problems of explaining the observed cosmic acceleration within GR, the implications of a consistent theory of modified gravity for cosmology are of immediate interest. Unfortunately, the original idea that a large vacuum energy contribution from the matter sector could be screened out simply by a non-zero graviton mass that weakens gravity at large scales turned out not to be realisable without fine tuning, as discussed in [46]. Nevertheless, one can take a less ambitious approach and assume that another, as of yet unknown, mechanism (such as supersymmetry, if it was realised at low energies) is at work and removes all vacuum energy. In this case, ordinary GR would predict a Universe without cosmic acceleration, while in modified gravity theories, ‘self-accelerating’ solutions certainly exist in various forms Massive gravity is a particularly interesting candidate for this purpose since its interaction potential breaks diffeomorphism invariance and therefore the spin-2 mass scale is expected to be protected from receiving large quantum corrections. The hope would therefore be that the interactions of the graviton could give rise to a small rate of cosmic acceleration, which, unlike the pure cosmological constant, can be regarded as being ‘technically natural’ in the sense of 'tHooft [198] (see e.g. [199] for a recent review on naturalness in the context of cosmology). Of course, it should be pointed out that although the guide of naturalness is philosophically appealing, in the end it may turn out that nature does not follow that principle. Nevertheless, like Occam’s razor, it is a useful tool in discriminating between alternative hypotheses.

Shortly after the construction of the ghost-free potential, it was discovered that the dRGT model of massive gravity (with flat or even with general reference metric) does not possess stable homogeneous and isotropic solutions [127]. For more work on cosmology in massive gravity with fixed $f_{\mu\nu}$, see [128–133, 200]. On the other hand, in bimetric theory with dynamical reference metric the desired class of solutions does exist [153–155]. The simplest examples are the proportional backgrounds considered in section 5.3.1 which, however, only solve the equations with proportional sources. In this case, the effective cosmological constant in (5.11) receives contributions not only through vacuum energy from the matter sector (as captured by $\beta_0$) but also from all terms in the interaction potential. Even in the absence of
vacuum energy (i.e., for $\beta_0 = 0$), cosmological backgrounds can be accelerating and the scale of acceleration is proportional to the technically natural mass scale $m$ (see e.g. [201]).

In the following we provide a brief summary of results in bimetric cosmology. For a more detailed discussion, we refer the reader to [202].

### 6.3.1. Homogeneous & isotropic backgrounds in bimetric theory.
We will outline the derivation of cosmological solutions in bimetric theory, following mainly [154]. For simplicity, the source $T_{\mu\nu}^f$ of the $f_{\mu\nu}$ sector is set to zero, which allows us to interpret $g_{\mu\nu}$ as the physical metric in the usual way, provided that it has standard matter couplings. In other analyses, the second source term is kept nonzero [158, 203, 204], which in principle can serve to mimic a dark matter component [149–152].

To keep the comparison of the bimetric solutions to those of GR as simple as possible, it is convenient to make the usual Friedmann–Lemaître–Robertson–Walker (FLRW) ansatz for the metric $g_{\mu\nu}$, which is coupled to matter$^{35}$

$$g_{\mu\nu}dx^{\mu}dx^{\nu} = -dt^2 + a(t)^2\left(\frac{1}{1 - kr^2}dr^2 + r^2d\Omega\right). \qquad (6.26)$$

It is now important to notice that, in order to arrive at this form, we have used time reparametrisations to remove one time dependent function in the most general homogeneous and isotropic ansatz. We have therefore fixed the entire diffeomorphism invariance and there is no gauge symmetry left to do the same to $f_{\mu\nu}$. As a result, the most general homogeneous and isotropic, bi-diagonal ansatz for the second metric reads

$$f_{\mu\nu}dx^{\mu}dx^{\nu} = -X(t)^2dt^2 + Y(t)^2\left(\frac{1}{1 - kr^2}dr^2 + r^2d\Omega\right), \qquad (6.27)$$

involving two time dependent functions $X(t)$ and $Y(t)$. Note that we have assumed the curvatures $k$ of the two metrics to be the same. This does not constitute a loss of generality because, starting with two different values for $k$, one can show that the Bianchi constraint forces them to be equal [205].

Furthermore, as in GR, we use the perfect fluid form $(T^g)^{\mu}_{\nu} = \text{diag}(-\rho, p, p, p)$ for the stress-energy tensor $T^g_{\mu\nu}$ of matter coupled to $g_{\mu\nu}$. In GR any source is automatically covariantly conserved as a consequence of the Bianchi identity whereas in bimetric theory the sources for $g_{\mu\nu}$ and $f_{\mu\nu}$ are not necessarily conserved. We therefore need to make the additional assumption of diffeomorphism invariance of the matter coupling which gives $\nabla^{\mu}T^g_{\mu\nu} = 0$. This in turn implies that the continuity equation for matter (see equation (C.4)) is also valid in the model under consideration.

Next, we consider the bimetric equations of motion in (6.1), with $T^f_{\mu\nu} = 0$ and with one index raised by the respective inverse metrics. In what follows we shall refer to their 00-components simply as the $g_{\mu\nu}$ and $f_{\mu\nu}$ equation. Plugging the ansätze of the previous section into (6.1a) we can make use of known GR results to straightforwardly obtain the $g_{\mu\nu}$ equation

$$\left(\frac{a}{a}\right)^2 + \frac{k}{a^2} - \frac{\mu^2}{3}\left[\beta_0 + 3\beta_1\frac{Y}{a} + 3\beta_2\left(\frac{Y}{a}\right)^2 + \beta_3\left(\frac{Y}{a}\right)^3\right] = \frac{\rho}{3m_g^2}, \qquad (6.28)$$

where $\mu^2 = m^4/m_g^2$. The first two terms are the same as in the ordinary Friedmann equation of GR, see (C.5). In order to determine how the source $\rho$ influences the cosmological evolution of the scale factor $a(t)$, we need to determine the function $Y(t)$ from the other equations. We start with the Bianchi constraint $\nabla^{\mu}V^g_{\mu\nu} = 0$, which, when evaluated on the

<footnote>$^{35}$ More general non-bidiagonal ansätze for the metrics have been studied in [153]. As we discussed in section 6.2.2 these are identical to GR backgrounds.</footnote>
homogeneous and isotropic ansatz, reduces to

$$ \frac{3\mu^2}{a} \left[ \beta_1 + 2\beta_2 \frac{Y}{a} + \beta_3 \left(\frac{Y}{a}\right)^2 \right] (\dot{Y} - \dot{a}X) = 0. \qquad (6.29) $$

In addition to these, we have to consider the $f_{\mu\nu}$ background equation whose form is slightly more complicated than (6.28) due to the presence of the additional function $X(t)$ in the ansatz for $f_{\mu\nu}$. Instead of presenting its complete form here, we first replace $X(t) = \dot{Y}/\dot{a}$, which is the dynamical solution to the Bianchi constraint above$^{36}$. Then the $f_{\mu\nu}$ equation becomes

$$ \left(\frac{\dot{a}}{Y}\right)^2 + \frac{k}{Y^2} - \frac{\mu^2}{3\alpha^2} \left[ \beta_1 \left(\frac{a}{Y}\right)^3 + 3\beta_2 \left(\frac{a}{Y}\right)^2 + 3\beta_3 \frac{a}{Y} + \beta_4 \right] = 0, \qquad (6.30) $$

where, as usual, $\alpha \equiv m_f/m_g$. Multiplying this equation by $Y^2/a^2$ and subtracting it from the $g_{\mu\nu}$ equation (6.28) yields an algebraic equation

$$ \frac{\beta_3}{3}\Upsilon^4 + \left(\beta_2 - \frac{\beta_4}{3\alpha^2}\right)\Upsilon^3 + \left(\beta_1 - \frac{\beta_3}{\alpha^2}\right)\Upsilon^2 + \left(\rho_* + \frac{\beta_0}{3} - \frac{\beta_2}{\alpha^2}\right)\Upsilon - \frac{\beta_1}{3\alpha^2} = 0, \qquad (6.31) $$

which we have expressed in terms of the dimensionless quantities

$$ \Upsilon = \frac{Y}{a}, \qquad \rho_* = \frac{\rho}{3m^4}. \qquad (6.32) $$

From the quartic polynomial equation (6.31) one can obtain a solution for $\Upsilon$ which can be inserted into the $g_{\mu\nu}$ equation (6.28). After that, the latter will only contain $a(t)$, $\dot{a}(t)$ and $\rho(t)$, and thus becomes a modified Friedmann equation.

### 6.3.2. Classification of solutions.
In general, the cosmological evolution equations derived above allow for several branches of solutions which have been categorised and studied in detail [148, 153–155, 201, 206–211]. Here we summarise the main results.

To start with, there are two options to satisfy the Bianchi constraint equation (6.29); either the first or the second bracket must vanish.

(I) The non-dynamical (algebraic) branch: If at least two out of $\beta_1$, $\beta_2$ and $\beta_3$ are non-vanishing, then a constant $b$ can be chosen such that $\beta_1 + 2b\beta_2 + b^2\beta_3 = 0$. In this case, an algebraic solution to the Bianchi constraint is $Y(t) = b a(t)$. Even though in this case the 00-components of the two metrics are not necessarily proportional, the contribution from the interaction potential in (6.28) reduces to a cosmological constant term. The $g_{\mu\nu}$ equation hence simply becomes the ordinary Friedmann equation (C.5) with a cosmological constant $\Lambda = m^2(\beta_0 + 3b\beta_1 + 3b^2\beta_2 + b^3\beta_3)$ and thus the background solutions for $g_{\mu\nu}$ are degenerate with those of GR. At the level of linear perturbations, several degrees of freedom appear without kinetic terms and are thus expected to be strongly coupled [129, 203, 212]. On top of that, this branch seems to give rise to a non-perturbative ghost [131] and a late-time instability for the tensor modes [107]. Due to these pathologies, most of the literature focusses on branch (II) solutions. Note also that branch (I) does not exist if two out of $\beta_1$, $\beta_2$ and $\beta_3$ are equal to zero.

(II) The dynamical branch: The alternative solution to the Bianchi constraint is $X(t) = \dot{Y}/\dot{a}$. In contrast to branch (I), this is an evolution equation and allows for much more general solutions that can be very different from GR. In particular, $\Upsilon$ is not constant on this branch but a function of the matter density $\rho$, determined by (6.31).

<footnote>$^{36}$ The other algebraic solution, obtained from setting the first bracket in (6.29) to zero, will be discussed below.</footnote>
The second branch further splits up into several subbranches, corresponding to different solutions of the quartic polynomial equation (6.31). These can be classified according to the evolution of the ratio of scale factors, $\Upsilon = Y/a$ [207].

(IIa) *Infinite branch*: This branch is characterised by an $\Upsilon$ that becomes infinitely large at early times and whose value decreases when moving forward in time.

(IIb) *Finite branch*: The solution is such that the ratio of scale factors $\Upsilon$ evolves towards a finite asymptotic value in the infinite past. This branch cannot be obtained in models with $\beta_1 = 0$ (assuming $\beta_0 = 0$) [207].

(IIc) *Exotic branches*: These are all other branches which generically lead to bouncing cosmologies or static universes in the asymptotic past or future [210]. Even though it is not impossible to make such scenarios compatible with observations, these branches correspond to rather unconventional cosmologies and have not been studied in great detail so far.

Using phenomenologically inspired requirements, one can immediately rule out a set of bimetric models based on their cosmological background evolutions. For instance, out of the models with only one non-vanishing $\beta_n$, only the $\beta_0$- and the $\beta_1$-model can give rise to a viable cosmology at the background level. Allowing for two non-vanishing $\beta_n$, one can generically discard models with $\beta_0 = \beta_1 = 0$. More details on the viability conditions and exclusion of bimetric models based on background cosmology can be found in [201, 207].

### 6.3.3. Self-acceleration & technical naturalness of the Hubble scale.
In view of the dark energy problem, it is of particular interest to investigate whether bimetric theory can give rise to an expansion history compatible with observations even when there is no vacuum energy arising from the matter sector. If a background solution can mimic the behaviour of vacuum energy, it is referred to as 'self-accelerating'. Clearly, it would be even more interesting if it was possible to include a large vacuum energy contribution and still obtain a small acceleration rate. This can indeed be achieved in bimetric theory by (fine-)tuning the $\beta_n$ parameters to cancel terms from the interaction potential against the vacuum contribution. However, there is no mechanism (or symmetry) protecting the tiny difference of these two contributions and thus a configuration with small effective cosmological constant cannot be regarded as technically natural in this setup. At this stage bimetric theory is therefore in no better position than GR to solve the old cosmological constant problem.

Nevertheless, it is still possible to view the spin-2 interaction potential as the source of dynamical dark energy, much in the same spirit as in quintessence models [213]. To this end one assumes that all vacuum contributions in the source $T_{\mu\nu}^g$ vanish and furthermore sets the bare cosmological constant term to zero by fixing $\beta_0 = 0$. A thorough comparison to observational data of the whole bimetric parameter space but with vanishing $\beta_0$ was performed in [201]. Self-accelerating solutions are found for many different combinations of nonzero $\beta_n$ parameters, and the fits to data are just as good as for the ordinary $\Lambda$CDM model of GR. As could have been anticipated, the best fit value for the mass of the massive spin-2 field is on the order of the Hubble scale $H_0$. This tiny mass scale creates the hierarchy problem in GR with an ordinary cosmological constant which is not protected against large quantum corrections by any symmetry. The situation is different in bimetric theory where in the zero-mass limit the full diffeomorphism symmetry that transforms the two metrics separately is restored. This symmetry protects the spin-2 mass scale against large corrections and renders
the small Hubble scale technically natural. For explicit calculations confirming this last statement, we refer the reader to [170, 214].

**6.3.4. Perturbations.** In order to further study the viability of bimetric cosmology, one derives the linear equations for perturbations around the above cosmological backgrounds. In contrast to what we encountered for the proportional backgrounds in section 5.3.1, it is not possible to diagonalise the fluctuations into spin-2 mass eigenstates around more general solutions. For the homogeneous and isotropic backgrounds, the best one can do is to perform a decomposition into tensor, vector and scalar modes and try to decouple their dynamics as much as possible. Rather general analyses of the spectrum were already performed in [212, 215] (and later in [216]), but since the resulting differential equations are not always easy to handle, developing a complete understanding of cosmological perturbation theory is still an ongoing process.

As already mentioned in section 6.3.2, branch (I) solutions are plagued by pathologies such as strong coupling behaviour and a non-perturbative ghost instability. Let us therefore focus on linear perturbations in the dynamical subbranches.

**Infinite branch (IIa):** While specific infinite branch solutions for models with $\beta_2 = \beta_3 = 0$ are free of scalar instabilities [208], in general the perturbations around this class of FLRW backgrounds suffer from ghosts in both the scalar and tensor sectors [209, 210, 217]. Due to these pathologies, branch IIa solutions do not correspond to consistent backgrounds and must be discarded on theoretical grounds. Nevertheless, for an interesting study of the integrated Sachs-Wolfe-effect in this branch which appeared before its problems were established, see [218].

**Finite branch (IIb):** Around this branch linear perturbations are generically ghost-free and well-behaved, except for a scalar gradient instability which occurs at early times [208, 212, 217, 219]. While, from a theoretical point of view, this growing scalar mode is not a consistency problem, it does invalidate the use of linear perturbation theory. For generic parameters, the instability sets in at recent times and, as a consequence, bimetric theory cannot be invoked to predict phenomena in the early Universe, at least not with standard perturbative techniques. However, it was demonstrated in [148] that, in the GR limit (see section 5.4.2), the scalar instability is pushed backwards to arbitrarily early times. Hence, for small enough $\alpha = m_f/m_g$, linear perturbation theory around the finite branches remains valid and the predictions for cosmology automatically resemble those made by GR. Another suggestion is that nonlinear effects related to the Vainshtein mechanism may render the instability irrelevant [220].

**Exotic branches (IIc):** As already mentioned, the background evolution on these branches is rather unusual and hence this type of solutions has not received much attention. The perturbations around these backgrounds generically seem to have pathologies [210].

To summarise, the only well-studied models that give rise to viable cosmological backgrounds and perturbations lie on the finite branch. In order to maintain the predictivity of cosmological perturbation theory, it is furthermore necessary to bring bimetric theory close to GR, either by requiring a sufficiently small value of $\alpha = m_f/m_g$ or by invoking the Vainshtein mechanism. The advantage of the bimetric setup is the occurrence of a technically natural dark energy scale set by the spin-2 mass. However, eventually all of this will be useful only if one finds an explanation for the absence of the large vacuum energy contribution coming from the matter sector.
## 7. Partial masslessness

A massive spin-2 field in de Sitter (dS) background exhibits several interesting features that are shared neither by lower-spin particles nor by spin-2 excitations in flat space. For a special value of the Fierz-Pauli mass in units of the background curvature, the linear theory possesses an additional gauge symmetry which removes one of the spin-2 helicity components. In this case, where the field loses one propagating degree of freedom, it is referred to as partially massless (PM). The particular mass value for which this situation occurs is called Higuchi bound and it divides the parameter region into unitary and non-unitary sub-sectors [81]. Below the Higuchi bound the helicity-zero mode of the spin-2 particle develops a ghost instability, whereas above the bound all helicity states are well-behaved.

The existence of a nonlinear theory that involves a massive spin-2 field enables us to address the question whether the concept of partial masslessness is restricted to the linear theory around de Sitter backgrounds or whether it may be extended to the nonlinear level. By definition, the demand on a nonlinear PM theory is the presence of an extra gauge symmetry, even away from de Sitter backgrounds.

### 7.1. PM spin-2 field on de Sitter background

Let us begin by reviewing the linear theory of a PM spin-2 particle in a de Sitter background, as first discussed in a sequence of papers by Deser et al [84–86, 194, 221–223]. Consider a de Sitter background metric $\bar{G}_{\mu\nu}$ whose curvature satisfies $R_{\mu\nu}(\bar{G}) = \Lambda\bar{G}_{\mu\nu}$ with positive cosmological constant $\Lambda$. The linearised action for a massive spin-2 fluctuation propagating on this background is

$$S_{\text{dS}} = \frac{1}{2} \int d^4 x \left[ h_{\mu\nu} \bar{\mathcal{E}}^{\mu\nu\rho\sigma} h_{\rho\sigma} - \Lambda \left( h^{\mu\nu} h_{\mu\nu} - \frac{1}{2} h^2 \right) + \frac{m_{\text{FP}}^2}{2} (h^{\mu\nu} h_{\mu\nu} - h^2) \right], \quad (7.1)$$

where the linearised Einstein operator $\bar{\mathcal{E}}^{\rho\sigma}_{\mu\nu}$ is the same as in (518). The corresponding equations of motion for $h_{\mu\nu}$ reads

$$\bar{\mathcal{E}}^{\rho\sigma}_{\mu\nu} h_{\rho\sigma} - \Lambda \left( h_{\mu\nu} - \frac{1}{2} \bar{G}_{\mu\nu} h \right) + \frac{m_{\text{FP}}^2}{2} (h_{\mu\nu} - \bar{G}_{\mu\nu} h) = 0. \quad (7.2)$$

For general values of the Fierz-Pauli mass, this equation possesses no gauge symmetries. As we showed in section 6.2.1, it gives rise to five constraint equations that reduce the number of propagating degrees of freedom to five, as appropriate for a massive spin-2 particle.

A remarkable feature of (7.2) without any analogue in flat space is that if the Fierz-Pauli mass saturates the Higuchi bound, $m_{\text{FP}}^2 = \frac{2}{3}\Lambda$, the equation is invariant under a new gauge symmetry. The corresponding linear gauge transformation reads

$$\Delta(h_{\mu\nu}) = \bar{\nabla}_{\mu} \partial_{\nu} \xi(x) + \frac{\Lambda}{3} \bar{G}_{\mu\nu} \xi(x), \quad (7.3)$$

in which $\xi(x)$ is a local gauge parameter.

Besides computing the variation of (7.2) under the gauge transformation explicitly, there is another straightforward way to see the existence of a gauge symmetry in the equation for the massive spin-2 field with mass at the Higuchi bound. First take the double divergence of (7.2) to arrive at

$$\frac{m_{\text{FP}}^2}{2} (\bar{G}^{\mu\rho} \bar{G}^{\nu\sigma} - \bar{G}^{\mu\nu} \bar{G}^{\rho\sigma}) \bar{\nabla}_{\nu} \bar{\nabla}_{\mu} h_{\rho\sigma} = 0. \quad (7.4)$$

The kinetic terms have dropped out after using the linearised Bianchi identity. Furthermore, the trace of (7.2) is given by