$$ (\bar{G}^{\mu\nu}\bar{G}^{\rho\sigma} - \bar{G}^{\mu\rho}\bar{G}^{\nu\sigma})\bar{\nabla}_{\nu}\bar{\nabla}_{\mu}h_{\rho\sigma} + (\Lambda - \frac{3}{2}m_{\text{FP}}^2)h = 0. \qquad (7.5) $$

If the mass is at the Higuchi bound, the terms without derivatives in this equation vanish identically, while the derivative terms are identical to those in (7.4). Hence we find that for $m_{\text{FP}}^2 = \frac{2}{3}\Lambda$ the traced equations of motion are identical to their double divergence or, in other words, the sum of the double divergence and the (correctly normalised) trace is identically zero. A gauge identity of this type even implies that the linearised action (7.1) is invariant under the corresponding gauge transformation. This can be seen by noticing that invariance of the action under (7.3) requires

$$ \Delta S_{\text{ds}} = \frac{\delta S_{\text{ds}}}{\delta h_{\mu\nu}} \Delta(h_{\mu\nu}) = \frac{\delta S_{\text{ds}}}{\delta h_{\mu\nu}} \left( \bar{\nabla}_{\mu} \partial_{\nu} \xi(x) + \frac{\Lambda}{3} \bar{G}_{\mu\nu} \xi(x) \right) = 0. \qquad (7.6) $$

Since $\frac{\delta S}{\delta h_{\mu\nu}}$ is nothing else than the equation of motion for $h_{\mu\nu}$, after integrating by parts, we find that $\Delta S$ indeed vanishes due to the gauge identity.

The conclusion is that, owing to the gauge redundancy, the spin-2 field described by (7.2) with mass at the Higuchi bound has only four dynamical degrees of freedom, one less than is the case for an ordinary massive spin-2 mode. The existence of the additional gauge symmetry around de Sitter background is in agreement with the representation theory of the de Sitter group SO (1, 4), which allows for ‘short’ representations of higher-spin fields. These representations contain less degrees of freedom than the massive ones and have been dubbed PM. For a discussion of this topic, consult e.g. [224].

An interesting feature of PM fields is that the enhanced symmetry could potentially improve the quantum behaviour of spin-2 theories. For instance, the gauge invariance protects the difference of cosmological constant and spin-2 mass against receiving large quantum corrections. Since the spin-2 mass scale itself is protected by diffeomorphism invariance of the massless theory, this could render a small value for the cosmological constant technically natural. Unfortunately, matter couplings of $h_{\mu\nu}$ destroy the gauge symmetry already at the linear level, unless the matter source is conformally invariant. An idea how to avoid this problem is presented in [225]. It is also worth mentioning that the linear PM theory possesses interesting properties similar to the electromagnetic duality in Maxwell’s theory [226, 227].

## 7.2. The search for a nonlinear PM theory

After the linear PM theory for spin-2 fields on de Sitter background had been discovered, it was soon attempted to construct higher-order interactions for the perturbation $h_{\mu\nu}$ that would leave the gauge symmetry intact. At cubic order in $h_{\mu\nu}$ around dS backgrounds such a construction turned out to be possible in $d = 4$ space-time dimensions [228–230]. On the other hand, in $d > 4$, it was found that gauge invariant cubic vertices can only exist if one also includes higher-derivative terms into the theory [228]. The construction of a fully nonlinear action with a PM gauge symmetry remains an open task. Many recent findings point towards the fact that such a theory cannot exist, in particular not as a theory involving nothing but a PM spin-2 field [134, 231–237]. We will summarise these arguments in section 7.4.

It is natural to expect the nonlinear PM theory (if it exists) to be found within the family of ghost-free nonlinear theories for massive spin-2 fields. A first investigation of this possibility was carried out in [238] where the authors aimed to identify a PM theory in nonlinear massive gravity with the reference metric taken to be a fixed de Sitter background. Their construction used Stückelberg fields in a generalised decoupling limit and showed that, for a certain choice of the $\beta_n$ parameters in the potential, a scalar degree of freedom is not
propagating. Here we will follow the slightly simpler strategy of [239] which derived the parameters for the PM candidate theory in bimetric theory$^{37}$. As it turns out, this particular bimetric model possesses several rather unexpected properties which could be of interest even if the no-go results against a nonlinear PM theory cannot be evaded. It is also possible that the nonlinear PM theory requires the input of additional degrees of freedom (e.g. higher-spin fields) and the framework of bimetric theory allows us to study merely the remnants of the enhanced symmetry present in the unknown extended setup.

All our considerations here will be in four space-time dimensions. It is possible to generalise some of the results to higher dimensions if the bimetric action is augmented by the Lanczos-Lovelock (LL) terms, which we shall come back to in section 8.1.

### 7.2.1. Identifying PM candidates.
Our first observation is that, in bimetric theory, the equation of motion (5.20b) for the massive fluctuation $\delta M_{\mu\nu}$ around proportional backgrounds is of the same form as (7.2) with $\Lambda = \bar{\Lambda}_g$. As a consequence, if the cosmological constant is assumed to be positive and if the Fierz-Pauli mass in (5.21) is on the Higuchi bound, a gauge symmetry of the form (7.3) is present in the linear theory around proportional backgrounds. In bimetric notation, the corresponding infinitesimal gauge transformation of the massive fluctuation reads

$$ \Delta(\delta M_{\mu\nu}) = \bar{\nabla}_{\mu} \partial_{\nu} \xi(x) + \frac{\bar{\Lambda}_g}{3} \bar{G}_{\mu\nu} \xi(x). \qquad (7.7) $$

On the other hand, the massless fluctuation transforms under the PM symmetry at most by a term that resembles a coordinate transformation, $\Delta(\delta G_{\mu\nu}) \sim \bar{\nabla}_{\mu} \partial_{\nu} \xi(x)$. This follows from the fact that its equation of motion (5.20a) does not have any extra gauge symmetry besides the usual linearised diffeomorphism invariance. Hence, we can always undo the PM transformation of $\delta G_{\mu\nu}$ by a coordinate transformation. Modulo GCTs, we can therefore demand $\Delta(\delta G_{\mu\nu}) = 0$.

The question is now under what conditions the above symmetry transformations for the fluctuations around dS background have a chance to be extendable to the nonlinear level. Since the nonlinear theory is formulated in terms of the variables $g_{\mu\nu}$ and $f_{\mu\nu}$, we first need to translate the transformations of mass eigenstates at the linear level into transformations of the fluctuations of the metrics. Assuming that we simultaneously perform a GCT to achieve $\Delta(\delta G_{\mu\nu}) = 0$, this translation can be obtained uniquely using the expressions for the massive and massless fluctuations given in (5.19). The result is

$$ \Delta(\delta g_{\mu\nu}) = -\alpha^2 \Delta(\delta f_{\mu\nu}) = \frac{-2(\alpha c)^2}{c(1+(\alpha c)^2)} \Delta(\delta M_{\mu\nu}). \qquad (7.8) $$

The crucial observation now is that, in a theory with a gauge symmetry at the nonlinear level, it must be possible to shift these transformations of the fluctuations $\delta g_{\mu\nu}$ and $\delta f_{\mu\nu}$ to the backgrounds $\bar{g}_{\mu\nu}$ and $\bar{f}_{\mu\nu}$. This follows simply from the fact that the split into background and fluctuations is not unique and one can always take out part of the fluctuation to redefine the background. Therefore we demand that the symmetry transformations (7.8) should also leave the background equations invariant. But since we are dealing with the proportional solutions, on which the background equations are reduced to an Einstein equation for $\bar{g}_{\mu\nu}$ along with $\Lambda_g = \Lambda_f$, see (5.12), we restrict ourselves to constant gauge transformations $\bar{\Delta}$ with

<footnote>$^{37}$ Note that the parameters of the PM model in massive gravity obtained in [238] are different from those first derived in [239] and later confirmed by [233, 235]. The reason for this mismatch can be traced back to a wrong choice of background, and correcting for this one can use the method of [238] to arrive at the parameters that are consistent with the other references.</footnote>
$\xi(x) = \xi_0 = \text{const.}$, such that $\Delta(\delta M_{\mu\nu}) = \frac{\xi_0 \bar{\Lambda}_g}{3} \bar{G}_{\mu\nu}$. In this way we ensure that the transformation does not take us away from the proportional backgrounds and avoid unnecessary complication. The restriction to constant gauge transformations is a strong simplification, but as we will see now, demanding the background equation to be invariant under these is constraining enough to identify the PM candidate theory. Namely, it should be evident from (7.8) that shifting the constant transformations to the backgrounds $\bar{g}_{\mu\nu}$ and $\bar{f}_{\mu\nu} = c^2 \bar{g}_{\mu\nu}$ results in a shift in $c^2$

$$ \bar{g}_{\mu\nu} \bar{\Delta}(c^2) = \bar{\Delta}(\bar{f}_{\mu\nu}) - c^2 \bar{\Delta}(\bar{g}_{\mu\nu}) = \bar{g}_{\mu\nu} \frac{2}{3} \Lambda_g c \xi_0. \qquad (7.9) $$

In general this cannot lead to an integrable<fnref n="38" /> symmetry of the background equations because $c$ is determined by $\Lambda_g = \Lambda_f$ and can therefore not be shifted. If this is the case, then a nonlinear PM symmetry cannot exist. The only possibility to avoid this immediate no-go statement is to demand that $c$ is not fixed by the background equation. Then, the equation $\Lambda_g = \Lambda_f$ that is automatically satisfied in this case (i.e. without fixing $c$) can be thought of as the gauge identity evaluated on the proportional backgrounds. Since the proportional ansatz for the two metrics partly fixes the gauge, only constant scalings are left as residual transformations.

### 7.2.2. PM candidate theory.
In four space-time dimensions, the explicit expressions (5.11) for the cosmological constants in terms of the $\beta_n$ parameters, the ratio of Planck masses $\alpha$ and the proportionality constant $c$ can be used to write the background condition $\Lambda_g = \Lambda_f$ as a polynomial equation for $c$

$$ \beta_1 + (3\beta_2 - \alpha^2\beta_0)c + 3(\beta_3 - \alpha^2\beta_1)c^2 + (\beta_4 - 3\alpha^2\beta_2)c^3 - \alpha^2\beta_3c^4 = 0. \qquad (7.10) $$

This equation clearly fixes $c$ unless all coefficients of different powers of $c$ vanish separately. A proportionality constant $c$ that is undetermined by the background equation therefore requires the following parameter choices in the bimetric interaction potential

$$ \beta_1 = \beta_3 = 0, \qquad \alpha^2\beta_0 = 3\beta_2 = \alpha^{-2}\beta_4. \qquad (7.11) $$

We will frequently refer to these values as the PM parameters. Note that our requirement on the $\beta_n$ parameters fixes all but one of them uniquely in terms of the others and $\alpha$. The one remaining parameter is of course degenerate with the scale $m^4$ in the interaction potential and sets the scale for the Fierz-Pauli mass and the cosmological constant. Moreover, it is easy to check using the expressions (5.21) and (517) for $\bar{m}_{\text{FP}}^2$ and $\bar{\Lambda}_g$, that the parameter choice (7.11) automatically puts the mass on the Higuchi bound. Therefore the theory of linear perturbations around the backgrounds at hand exhibits the usual PM gauge symmetry. It is worth emphasising that we did not demand this in any way; it followed from an independent requirement on the background equations<fnref n="39" />. Starting from the proportional backgrounds the finite form of the scaling symmetry is given, for any constant $a$, by [239]

<footnote n="38">³⁸ Integrability of the symmetry transformation means that it can be performed more than once and still leave the equations invariant. This is a natural requirement on any gauge symmetry. Without the integrability condition it is sufficient to set the mass to the Higuchi bound in order to have an invariant background.</footnote>

<footnote n="39">³⁹ On the other hand, note that if we assume that there is a unique PM theory, then requiring $\bar{m}_{\text{FP}}^2 = \frac{2}{3}\bar{\Lambda}_g$ is sufficient to determine its parameters. Namely, in this case, the symmetry (5.6) of the interaction potential enforces $\alpha^{4-n}\beta_n = \alpha^n\beta_{4-n}$ because otherwise the theory obtained from replacing $\alpha^{4-n}\beta_n \to \alpha^n\beta_{4-n}$ would also be PM, contradicting the uniqueness requirement. It is easy to see that this constraint on the $\beta_n$ parameters together with the Higuchi bound condition already implies (7.11).</footnote>
$$c \to c + a, \qquad \bar{g}_{\mu\nu} \to \frac{1 + (\alpha c)^2}{1 + \alpha^2 (a + c)^2} \bar{g}_{\mu\nu}. \qquad (7.12)$$

By restricting to constant gauge transformations, clearly, we are not dealing with the full PM symmetry at the nonlinear level. Note, however, that since the set of constant transformations is a subset of the full gauge group, the theory we obtain by requiring invariance under this subset must contain the theory with the full gauge group (if it exists). Moreover, since in four-dimensions the restriction to constant gauge transformations is sufficient to uniquely determine the $\beta_n$ parameters, we conclude that the resulting theory is already the unique candidate for having the full PM symmetry.

The scale invariance of the equations of motion for proportional backgrounds is not the only interesting property of the theory specified by the PM parameters. The PM candidate exhibits additional astonishing features that further support the existence of a gauge symmetry at the nonlinear level. For instance, consider again the homogeneous and isotropic solutions which we presented in section 6.3.1. We saw that, after solving the Bianchi constraint, it was possible to arrive at (6.31), an algebraic equation for the ratio of the two scale factors, $Y(t)/a(t)$. Inserting the PM parameters (7.11) into this equation in vacuum, we find that it becomes an identity. Clearly, it is the analogue of the equation $\Lambda_g = \Lambda_f$ for proportional backgrounds. But instead of a constant, the equations evaluated on the cosmological ansatz now leave a time-dependent function undetermined. This shows that the cosmological evolution equations of PM bimetric theory are invariant under a symmetry that is local in time.

Since, in four-dimensions, the proportional backgrounds correspond to maximally symmetric spacetimes with ten independent Killing vectors and the homogeneous and isotropic solutions still possess six isometries, one might speculate that the presence of a symmetry on these solutions could somehow be related to the amount of symmetry of the underlying geometry. To obtain more general results, it is therefore important to investigate the structure of our PM candidate theory beyond the proportional and cosmological backgrounds. We shall come back to this point in the next subsection where we show that the equations in the PM theory are Weyl invariant to lowest order in a derivative expansion.

**Massive gravity limit.** Let us comment briefly on the massive gravity limit. The PM parameters (7.11) in four space-time dimensions provide an example for a theory in which the coefficients in the interaction potential depend on the ratio $\alpha = m_f/m_g$ that is taken to infinity in the massive gravity limit which we discussed in section 5.4.3. In this case, we have to be more careful when taking the limit because the $\beta_n$ parameters will scale as well. The conditions (7.11) fix the relative scale among them but their absolute scale may still be chosen freely. Suppose that $\beta_2$ does not scale when the limit $\alpha = m_f/m_g \to \infty$ is taken. Then $\beta_0 \to 0$ and $\beta_4 \to \infty$, whereas $\Lambda = \beta_4 m^4/m_f^2 = \text{const}$. In the massive gravity limit, the equations of motion for the PM theory thus reduce to

$$G_{\mu\nu}(g) + \frac{\Lambda}{3} g_{\mu\lambda} (Y^{(2)})^{\lambda}_{\nu}(S) = 0, \qquad G_{\mu\nu}(f) + \Lambda f_{\mu\nu} = 0. \qquad (7.13)$$

These are exactly the equations singled out in [233, 235] which investigated the possibility of realising nonlinear partial masslessness in massive gravity with fixed reference metric. The methods invoked in those references are completely different from ours and could provide independent support for the existence of a nonlinear PM symmetry. On the other hand, as we already mentioned previously, [233, 235] also provided evidence for the absence of an additional gauge symmetry present in (7.13). These results may extend to the bimetric case
but such a generalisation is not obvious due to the singular nature of the massive gravity limit [147].

### 7.3. Connection to conformal gravity

As shown in [240], the PM candidate model is closely related to the well-known theory of conformal gravity whose action is invariant under Weyl transformations of the metric. In the following we will briefly review conformal gravity and then summarise the steps that establish its connection to the PM bimetric model.

#### 7.3.1. Review of conformal gravity.
In four space-time dimensions, there exists a particular higher-derivative action with an additional gauge symmetry [241]

$$S_{\text{CG}} = -\sigma \int d^4 x \sqrt{g} \left( R^{\mu\nu} R_{\mu\nu} - \frac{1}{3} R^2 \right), \qquad (7.14)$$

with dimensionless coefficient $\sigma$. This action is invariant under Weyl transformations of the metric,

$$g_{\mu\nu} \mapsto \xi^2(x) g_{\mu\nu}, \qquad (7.15)$$

where $\xi(x)$ is a local gauge parameter. The cosmological constant and the Einstein-Hilbert term are not invariant under the transformation (7.15); hence they do not appear in the above conformal gravity action.

The equations of motions that follow from variation of the conformal gravity action (7.14) imply the vanishing of the Bach tensor for $g_{\mu\nu}$ [242]

$$B_{\mu\nu} \equiv -\nabla^2 P_{\mu\nu} + \nabla^\rho \nabla_{(\mu} P_{\nu)\rho} + W_{\rho\mu\nu\sigma} P^{\rho\sigma} = 0. \qquad (7.16)$$

Here, we have given the definition of $B_{\mu\nu}$ in terms of the Schouten tensor

$$P_{\mu\nu} \equiv R_{\mu\nu} - \frac{1}{6} g_{\mu\nu} R, \qquad (7.17)$$

as well as the Weyl tensor$^{40}$

$$W_{\rho\mu\nu\sigma} \equiv R_{\rho\mu\nu\sigma} + g_{\mu\nu} R_{\sigma\rho} - g_{\rho\nu} R_{\sigma\mu} + \frac{1}{3} g_{\rho\nu} g_{\sigma\mu} R. \qquad (7.18)$$

The Bach equation (7.16) is also invariant under Weyl transformations in four space-time dimensions.

The conformal gravity action is closely related to the linear theory for a PM spin-2 field discussed in section 7.1. In order to see this, one can introduce an auxiliary field $\varphi_{\mu\nu}$ and a parameter $\Lambda$ to write down an equivalent action [243],

$$S_{\text{CG}}^{(2)} = 4\sigma\Lambda \int d^4 x \sqrt{g} \left[ \frac{1}{6}(R - 2\Lambda) + \varphi^{\mu\nu} G_{\mu\nu} + \Lambda\varphi + \Lambda(\varphi^{\mu\nu}\varphi_{\mu\nu} - \varphi^2) \right]. \qquad (7.19)$$

Linearisation around a constant curvature background and diagonalisation into mass eigenstates shows that, for $\sigma > 0$, the theory propagates a healthy massless and a ghost-like massive spin-2 particle. Moreover, the Fierz-Pauli mass of the massive fluctuation has the value corresponding to the Higuchi bound, $m_{\text{FP}}^2 = \frac{2}{3}\Lambda$, which implies the presence of the PM gauge symmetry (7.3) for $\varphi_{\mu\nu}$. The massless fluctuation does not transform under this symmetry. In fact, it is a simple exercise to check that, in addition to the obvious

<footnote>$^{40}$ The conformal gravity action (7.14) can also be expressed, modulo the Euler invariant, as the square of the Weyl tensor, $L_{\text{CG}} \propto \sqrt{g} W_{\rho\mu\nu\sigma} W^{\rho\mu\nu\sigma}$.</footnote>
diffeomorphism invariance, the full nonlinear auxiliary action (7.19) is invariant under the following linear gauge transformations [231]

$$ \delta g_{\mu\nu} = 2\xi(x)g_{\mu\nu}, \quad \delta\varphi_{\mu\nu} = (\nabla_{\mu}\partial_{\nu} + \frac{\Lambda}{3}g_{\mu\nu})\xi(x). \qquad (7.20) $$

Note that the nonlinear field $g_{\mu\nu}$ transforms under the conformal part of the PM symmetry because its fluctuation does not correspond to the massless mode but is a linear superposition of mass eigenstates. Its transformation is of course nothing but the linear version of the Weyl transformation (7.15).

Conformal gravity, or its equivalent form (7.19), therefore describes only six propagating degrees of freedom instead of seven that would correspond to a massless and a massive spin-2 field. For more detailed discussions of its spectrum we refer the reader to [244–247]. Like any other Weyl invariant theory, $S_{CG}$ does not contain any dimensionful couplings and therefore avoids the non-renormalisibility problem of GR. Being a renormalisable field theory, conformal gravity has been suggested as a quantum theory for gravity. It has also been shown to allow for viable cosmological solutions and even fit Galaxy rotation curves, providing a possible solution for at least part of the dark matter problem [248].

Unfortunately, all of its features remain irrelevant unless a cure for the ghost problem in conformal gravity is found. Suggestions for altering the theory in order to make it healthy include, for instance, a modification of quantum mechanics [249] and the addition of specific boundary conditions [247], but none of these have been sufficiently convincing for the theory to be accepted as a consistent alternative to GR. In order to avoid the inconsistencies that plague theories with a finite number of higher derivatives, another possible solution is to complete the equations of motion (7.16) of conformal gravity with an infinite number of derivatives for which Ostrodgradski’s theorem does not hold. Higher-derivative terms will, of course, break the conformal symmetry because they necessarily enter with a suppressing mass scale. On the other hand, one could imagine that the symmetry transformation needs to be completed with an infinite number of higher-derivative terms as well, in order to be a gauge symmetry of the full theory. In that case, both the equations of motion and the symmetry transformation could be thought of as a perturbative expansion in derivatives. Order by order, the equations would be invariant under the gauge symmetry, starting with the Bach tensor and its Weyl invariance at lowest order. Although this idea sounds promising, without any further input it is difficult to guess the form of the higher-derivative corrections to the Bach equation that could give rise to such a gauge symmetry. It would therefore be helpful to have a guideline telling us how to construct these terms This is where the PM candidate of bimetric theory comes into play.

### 7.3.2. Perturbative solution to the $g_{\mu\nu}$ equation.
We review here the results of [240]. Our aim is to combine the bimetric equations of motion to eliminate one of the metrics and derive an effective equation involving only the other. To this end, we will solve the $g_{\mu\nu}$ equation algebraically for the square-root matrix $S^{\mu\nu}$ as a perturbation series in curvatures of $g_{\mu\nu}$. From this we deduce a perturbative solution for $f_{\mu\nu}$ that can be used to eliminate $f_{\mu\nu}$ from its own equation of motion, resulting in a perturbative equation for $g_{\mu\nu}$. We will derive the lowest orders of this equation for general $\beta_n$ parameters and then see that it exhibits remarkable features when we restrict to the PM parameters (7.11) in the subsequent subsection. Of course, we could switch the roles of the metrics and in a similar manner derive an effective equation for $f_{\mu\nu}$.
It will prove convenient to first rewrite the Einstein tensor in terms of the Schouten tensor defined in (7.17) and raise one index with $g^{\mu\nu}$. This gives$^{41}$

$$ \mathrm{Tr} P \delta_{\nu}^{\mu} - P^{\mu}_{\nu} = \mu^2 \sum_{n=0}^{3} (-1)^n \beta_n (Y^{(n)})^{\mu}_{\nu}(S), \quad \mu^2 \equiv \frac{m^4}{m_g^2}. \qquad (7.21) $$

Then we make the following perturbative ansatz for the square-root matrix $S = \sqrt{g^{-1}f}$

$$ S^{\mu}_{\nu} = a \delta^{\mu}_{\nu} + \frac{1}{\mu^2} (b_1 P^{\mu}_{\nu} + b_2 \mathrm{Tr} P \delta^{\mu}_{\nu}) + \frac{1}{\mu^4} (c_1 (P^2)^{\mu}_{\nu} + c_2 \mathrm{Tr} P P^{\mu}_{\nu} + c_3 \mathrm{Tr} P^2 \delta^{\mu}_{\nu} + c_4 (\mathrm{Tr} P)^2 \delta^{\mu}_{\nu}) + \mathcal{O}(\frac{P^3}{m^6}), \qquad (7.22) $$

with arbitrary complex coefficients $a, b_i, c_i, ...$. Allowing the parameters to assume complex values is reasonable as long as the coefficients that will finally appear in the effective equation for $g_{\mu\nu}$ remain real.

In the model where only $\beta_0, \beta_1$ and $\beta_4$ are non-vanishing, the solution takes the very simple closed form

$$ S^{\mu}_{\nu} = -\frac{\beta_0}{3\beta_1}\delta^{\mu}_{\nu} + \frac{1}{\beta_1\mu^2}P^{\mu}_{\nu}(g), \quad \text{if } \beta_2 = \beta_3 = 0. \qquad (7.23) $$

For more general parameters, the expansion does not terminate and we can only determine the coefficients in the ansatz order by order in curvatures. As a side-remark: the existence of this exact solution in the $\beta_1$ model is rather remarkable. For example, linearising the equations of motion in this model we can use the above relation to fully remove any occurrence of $f_{\mu\nu}$ in the linearised equations. In particular, in the massive gravity limit this means that it is possible to obtain equations for a massive spin-2 field propagating on any background without any reference to a second metric [104, 105].

In order to simplify the expressions in the following, we introduce a new set of linear combinations of the $\beta_n$ parameters

$$ s_n \equiv \sum_{k=n}^{3} \binom{3-n}{k-n} a^k \beta_k, \qquad (7.24) $$

Note that on proportional backgrounds, $f_{\mu\nu} = c^2 g_{\mu\nu}$ and $s_0$ is proportional to the cosmological constant $\Lambda_g$ defined in (5.11) if in that expression one replaces $c$ by $a$.

The lowest orders in the solution for $S$ are determined to be of the following form

$$ S^{\mu}_{\nu} = a \delta^{\mu}_{\nu} + \frac{a}{s_1 \mu^2} P^{\mu}_{\nu} + \frac{a s_2}{s_1^3 \mu^4} \left[ (P^2)^{\mu}_{\nu} - \mathrm{Tr} P P^{\mu}_{\nu} + \frac{1}{3} e_2(P) \delta^{\mu}_{\nu} \right] + \mathcal{O}\left(\frac{P^3}{\mu^6}\right). \qquad (7.25) $$

We can then use $f_{\mu\nu} = g_{\mu\rho} (S^2)^{\rho}_{\nu}$ to arrive at

$$ f_{\mu\nu} = a^2 g_{\mu\nu} + \frac{2a^2}{s_1 \mu^2} P_{\mu\nu} + \frac{a^2(s_1 + 2s_2)}{s_1^3 \mu^4} (P^2)^{\mu}_{\nu} + \frac{2a^2 s_2}{s_1^3 \mu^4} \left[ \frac{1}{3} e_2(P) \delta^{\mu}_{\nu} - \mathrm{Tr} P P^{\mu}_{\nu} \right] + \mathcal{O}\left(\frac{P^3}{\mu^6}\right). \qquad (7.26) $$

This is the most general covariant solution for $f_{\mu\nu}$ obtained from the $g_{\mu\nu}$ equation. An immediate consequence of equation (7.26) is that if the solution for $g_{\mu\nu}$ has constant curvature, $P_{\mu\nu} \propto g_{\mu\nu}$, then $f_{\mu\nu} \propto g_{\mu\nu}$, i.e. the two metrics are proportional to each other. But, as we already mentioned in section 6.1, there exist solutions for which both metrics have

<footnote>$^{41}$ Mainly for the sake of notational simplicity, we have set the source for $g_{\mu\nu}$ to zero. This simplification constitutes no loss of generality because in the final results of our computation it can always be reinstated by making the replacement, $P^{\mu}_{\nu} \rightarrow P^{\mu}_{\nu} - m_g^{-2} (T g^{\mu\nu} - \frac{1}{3} \mathrm{Tr} T g^{\mu\nu})$.</footnote>
constant curvature while not being proportional to each other. This class of (non-covariant) solutions is therefore not captured by our ansatz for S in (7.22).

Next we use the solution (7.26) for $f_{\mu\nu}$ to eliminate it in its own equation of motion. This will lead to a set of effective equations for $g_{\mu\nu}$ containing an infinite number of derivatives, due to the presence of the inverse $f^{\mu\nu}$ in these equations. Inserting (7.26) into the Einstein tensor for $f_{\mu\nu}$ we find

$$ G_{\mu\nu}(f) = G_{\mu\nu}(g) - \frac{1}{s_1\mu^2}(\nabla^2 P_{\mu\nu} + \nabla_{\mu}\nabla_{\nu}P - \nabla^{\rho}\nabla_{\mu}P_{\rho\nu} - \nabla^{\rho}\nabla_{\nu}P_{\rho\mu} + 3PP_{\mu\nu} - g_{\mu\nu}[P^{\alpha\beta}P_{\alpha\beta} + \frac{1}{2}P^2]) + \mathcal{O}(\frac{P^3}{\mu^4}). \quad (7.27) $$

The contributions from the interaction potential to the $f_{\mu\nu}$ equation evaluated on (7.26) become

$$ \frac{\mu^2}{\alpha^2} \tilde{V}_{\mu\nu} = \frac{\mu^2 \Omega}{a^2 \alpha^2} g_{\mu\nu} + \frac{1}{a^2 \alpha^2} G_{\mu\nu} + \frac{2\Omega}{a^2 \alpha^2 s_1} P_{\mu\nu} + \frac{1}{a^2 \alpha^2 s_1^3 \mu^2} [c_1 P_{\mu}^{\rho} P_{\rho\nu} + c_2 P P_{\mu\nu} + \frac{1}{6} g_{\mu\nu} (c_3 P^{\alpha\beta} P_{\alpha\beta} - c_2 P^2)] + \mathcal{O}(\frac{P^3}{\mu^4}), \quad (7.28) $$

in which the expansion coefficients are given by

$$ c_1 = 2s_1^2 + \Omega(s_1 + 2s_2), \quad c_2 = -3s_1^2 - 2s_2\Omega, \quad c_3 = 3s_1^2 - 2s_2\Omega, \quad (7.29) $$

where we have defined

$$ \Omega = a\beta_1 + 3a^2\beta_2 + 3a^3\beta_3 + a^4\beta_4. \quad (7.30) $$

Note that this would be proportional to $\Lambda_f$ in (5.11) if $f_{\mu\nu} = a^2 g_{\mu\nu}$. Combining the kinetic and potential terms, we can write the entire $f_{\mu\nu}$ equation of motion as a higher-derivative equation for $g_{\mu\nu}$

$$ \frac{\Omega}{a^2 \alpha^2} g_{\mu\nu} + \frac{1}{\mu^2} \left[ 1 + \frac{1}{a^2 \alpha^2} \right] G_{\mu\nu} + \frac{2\Omega}{a^2 \alpha^2 s_1 \mu^2} P_{\mu\nu} + \frac{1}{\mu^4 s_1} B_{\mu\nu} \\ + \frac{\Omega}{a^2 \alpha^2 s_1^3 \mu^4} \left[ (s_1 + 2s_2) P_{\mu}^{\rho} P_{\rho\nu} - 2s_2 P P_{\mu\nu} - \frac{s_2}{3} g_{\mu\nu} (P_{\rho\sigma} P^{\rho\sigma} - P^2) \right] \\ - \frac{1}{s_1 \mu^4} \left( 1 + \frac{1}{\alpha^2 a^2} \right) \left[ 3 P P_{\mu\nu} - 2 P_{\mu}^{\rho} P_{\rho\nu} - \frac{1}{2} g_{\mu\nu} (P^2 - P^{\alpha\beta} P_{\alpha\beta}) \right] \\ + \mathcal{O}\left(\frac{P^3}{\mu^6}\right) = 0, \quad (7.31) $$

where we have collected some of the terms with four derivatives into the Bach tensor defined in (7.16). We have thus arrived at an effective equation involving only $g_{\mu\nu}$ in which the terms with more than four derivatives can be computed order by order following the same procedure as outlined above. The terms in the first line of (7.31) are the cosmological constant, the Einstein tensor and a correction proportional to the Schouten tensor. The equation hence reduces to GR in the small curvature limit only if $(\Omega/a^2\alpha^2)(1/s_1\mu^2) \to 0$. Note that the first of these brackets directly sets the size of the observed cosmological constant. The phenomenological relevance of this equation for weak gravitational fields has not been studied so far, but doing so would require adding a source in accordance with our remark in footnote 38.

Before coming to the PM case, let us make one more remark: plugging the perturbative solution (7.26) back into the bimetric action results in an effective action of the form
$$S_{\text{HD}} = m_g^2 \int d^4 x \sqrt{g} \left[ c_{\Lambda} + c_R R(g) - \frac{c_{RR}}{m^2} \left( R^{\mu\nu} R_{\mu\nu} - \frac{1}{3} R^2 \right) \right] + \mathcal{O}\left( \frac{p^3}{m^6} \right), \qquad (7.32)$$

where $c_Λ$, $c_R$ and $c_{RR}$ are functions of the bimetric parameters. We recognise the conformal gravity combination at fourth order in derivatives and notice that the above action, at quadratic order in curvatures, represents the generalisation of three-dimensional New Massive Gravity [250] to four-dimensions which is also known to propagate seven degrees of freedom [241, 251]. However, the action (7.32) is not equivalent to the original bimetric action since we have used the equation for $g_{\mu\nu}$ instead of the one for $f_{\mu\nu}$ to obtain a solution for $f_{\mu\nu}$. The equations derived from the above action differ from the bimetric equations by the product of a differential operator. Only by restricting to solutions where the zero modes of this operator are absent will the equations give the same solutions. The additional terms that would arise using the correct procedure, i.e. integrating out $f_{\mu\nu}$ by its own equations of motion, are nonlocal because a derivative operator needs to be inverted in order to solve the $f_{\mu\nu}$ equation for $f_{\mu\nu}$. This is discussed in detail in [240, 252]. In the following, we will not work with any effective action but restrict ourselves to the equations of motion where this ambiguity does not arise.

### 7.3.3. PM case.
Let us determine the set of parameter values for which the higher-derivative equation (7.31) obtained from bimetric theory reduces to the Bach equation (7.16) at lowest order in derivatives. This means that we require

$$\Omega = 0, \qquad a^2 = -\alpha^{-2}. \qquad (7.33)$$

In addition to this, satisfying the $g_{\mu\nu}$ equation at lowest order in curvatures requires that the combination $s_0$ defined in (7.24) vanishes. Using $a^2 = -\alpha^{-2}$ in the expressions for $\Omega$ and $s_0$, it is easy to see that these requirements combine into two complex equations

$$\begin{aligned} & \beta_0 + \frac{3i}{\alpha}\beta_1 - \frac{3}{\alpha^2}\beta_2 - \frac{i}{\alpha^3}\beta_3 = 0, \\ & \beta_0 + \frac{4i}{\alpha}\beta_1 - \frac{6}{\alpha^2}\beta_2 - \frac{4i}{\alpha^3}\beta_3 + \frac{1}{\alpha^4}\beta_4 = 0, \end{aligned} \qquad (7.34)$$

which need to be solved for the $\beta_n$ parameters. The real and imaginary parts must vanish separately and the unique solution is

$$\beta_1 = \beta_3 = 0, \qquad \alpha^2\beta_0 = 3\beta_2 = \alpha^{-2}\beta_4. \qquad (7.35)$$

Remarkably, these values precisely corresponds to the PM parameter choice in (7.11). We conclude that the PM candidate is the unique member of the family of bimetric theories with the feature that to lowest order in a derivative expansion the effective equation for $g_{\mu\nu}$ is identical to the Weyl invariant equation of conformal gravity

$$B_{\mu\nu} + \mathcal{O}\left(\frac{p^3}{m^6}\right) = 0. \qquad (7.36)$$

Moreover, it is straightforward to show that the method that we used here to arrive at an effective equation for $g_{\mu\nu}$ can yield a similar result for $f_{\mu\nu}$: the effective equation for $f_{\mu\nu}$, obtained from solving the $f_{\mu\nu}$ equation for $g_{\mu\nu}$ and plugging the solution into the $g_{\mu\nu}$ equation, also sets the Bach tensor for $f_{\mu\nu}$ to zero at lowest order in derivatives. This can directly be deduced from the symmetry property (5.7) of the equations of motion and the invariance of the PM parameter choice under $\alpha^{4-n}\beta_n \rightarrow \alpha^n\beta_{4-n}$.

In this way, the nonlinear PM bimetric theory proposes a ghost-free completion of conformal gravity. The presence of a gauge symmetry to lowest order in a derivative expansion can be regarded as further support of the existence of an additional gauge
symmetry in the theory. In particular, note that the above analysis is not based on de Sitter (nor FLRW) background.

If existent, the PM symmetry could be viewed as the generalisation of the gauge invariance of conformal gravity to higher orders in derivatives. In principle, it is possible to explicitly demonstrate the invariance of the equations under symmetry transformations that are extended to higher orders in derivatives. Using the perturbative solution to the $g_{\mu\nu}$ equation together with the analogous expression obtained from the $f_{\mu\nu}$ equation, the gauge invariance has been shown to exist up to sixth order in derivatives [253].

### 7.4. Arguments against a nonlinear PM theory

As mentioned before, the literature contains a large variety of no-go theorems that forbid the existence of a nonlinear theory for PM spin-2 fields. Most of these counter arguments refer to an action involving no other degrees of freedom besides the PM field but some of them may eventually also apply to bimetric theory. Let us summarise them briefly:

* In the massive gravity limit of the PM bimetric model, first obstructions were already encountered by the authors of [235] and it was shown in [232, 233] that the equations of motion do not satisfy the nonlinear Bianchi identity which is expected in a PM theory. More explicitly, the authors showed that the nonlinear covariant constraint which removes the Boulware–Deser ghost does not identically vanish. This result was argued to extend to bimetric theory [234] where, however, the situation is less clear since the nonlinear version of the constraint is not known. In fact, a covariant constraint does not even seem to exist for general backgrounds [108].

* The linear spectrum of conformal gravity always contains a ghost [231]. Although the full nonlinear bimetric theory is ghost-free, it is possible that the appearance of the Weyl invariant Bach tensor (7.36) in our perturbative approach is intimately related to the fact that we are expanding in small curvatures, i.e. around flat solutions which suffer from the same pathology as conformal gravity. Flat backgrounds in the PM theory require choosing $c^2 = -1$, which renders the kinetic term of the massive perturbation ghost-like [253]. Hence, the Weyl invariance at lowest order in the equations of motion (7.36) seems to be obtainable only at the cost of giving up unitarity.

* The authors of [134] argued that, in order to give rise to a nonlinear PM symmetry, the scalar mode of the massive field needs to disappear entirely from a 'decoupling limit' of bimetric theory. It was then demonstrated that, even though the pure scalar interactions are indeed absent, the mode reappears in interactions with the vectors$^{42}$. The fact that the vectors vanish on maximally symmetric backgrounds explains why one sees a PM gauge symmetry to linear order around dS space but not beyond.

* Imposing a closure condition on the PM transformations, the authors of [237] showed that no nonlinear Lagrangian involving at most two derivatives on the fields can realise the required symmetry for a single spin-2 field in a gravitational background. The setup is more general than ghost-free massive gravity but it still does not include the bimetric case.

<footnote n="42">42 From our point of view, it is not obvious if this analysis is actually able to rule out the gauge symmetry because, firstly, the backgrounds that [134] assumed for the metrics do not solve the equations of motion and, secondly, the constraint that removes the ghost is not imposed. As a consequence of the latter, the vector-scalar interactions involve higher-derivative terms of the scalar mode that are known to disappear after accounting for the constraint.</footnote>
* More group theoretical evidence against a nonlinear PM symmetry for spin-2 fields coupled to gravity was provided in [236] where it was shown that no unitary theory exists and the unique non-unitary example is conformal gravity. However, for technical reasons which are discussed in [253], these arguments do not directly carry over to bimetric theory.

* In [254] the authors provide arguments against the existence of a PM theory with a non-abelian Yang-Mills like extension. Again, these arguments cannot directly be applied to the bimetric setup.

If the above results eventually turn out to also be extendable to bimetric theory, this would imply that the PM symmetry can indeed not be realised as a nonlinear theory involving only spin-2 fields. But even in this case, the bimetric candidate model seems to possess interesting properties that could pave the way towards understanding partial masslessness from a background-independent point of view. A promising future direction could be the combination of bimetric theory with higher-spin degrees of freedom. Finally, let us point out that finding an additional scalar symmetry is not merely an interesting exercise but could indeed prove to be very useful from several perspectives. First of all, it would guarantee that the helicity-zero mode is absent even nonlinearly. This will have an enormous effect on the phenomenology since the scalar mode is usually responsible for various instability issues and causes the most tension with observations. Secondly it would provide an argument for why a small cosmological constant is technically natural since its value is tied to the mass of the spin-2 field via a symmetry, while a small mass is itself technically natural since its vanishing restores full diffeomorphism invariance. Thirdly, the absence of the scalar mode may lead to improved quantum behaviour of the theory. Indeed, since the symmetry if it exists is given by Weyl scalings to lowest order in a derivative expansion, the theory may have much better renormalisation properties than GR.

# 8. Extensions of bimetric theory

Extensions of massive gravity have been proposed mainly by introducing new scalar degrees of freedom and in the context of cosmology. For instance, a scenario in which the spin-2 mass is obtained through a scalar condensate has been suggested in [127] and a quasi-dilaton extension of the massive gravity action was constructed in [255, 256]. Similarly $f(R)$ extensions have been studied in e.g. [138–142] and other scalar fields with non-minimal coupling terms have also been added to the bimetric action, as in [257].

In the following we shall focus on more direct generalisations which preserves the structure of the bimetric theory itself without adding additional degrees of freedom, except the most natural generalisation to include interactions of multiple massive spin-2 fields$^{43}$. First, we discuss its generalisation to higher-dimensions, where new kinetic terms satisfy the requirement of classical consistency. Then we shall review interactions of multiple spin-2 fields and, in the same context, also present the vierbein formulation of bimetric theory.

## 8.1. Higher derivative LL extension

According to the constructive consistency proof in section 3.4, the interaction potential of the Hassan-Rosen bimetric theory includes all ghost-free non-derivative terms. Therefore, the only option to obtain an extended version of the theory is to add more derivative terms to the

<footnote>$^{43}$ Interactions between multiple massless spin-2 fields are forbidden on quite general grounds [258].</footnote>
action. The literature contains several no-go theorems on generalising the kinetic structure of the spin-2 fields in four-dimensions to anything beyond the Einstein-Hilbert term<fnref n="44" /> [135, 136, 259]. On the other hand, in dimensions greater than four, the LL invariants, which either vanish or are topological in d = 4, are expected not to reintroduce the Boulware-Deser ghost [263–267]. These are defined as totally antisymmetric contractions of Riemann tensors R$^{\alpha\beta}$$_{\mu\nu}$

$$
\mathcal{L}_{(n)} = \frac{1}{2^n n!} \delta^{\mu_1}_{[\alpha_1} \delta^{\nu_1}_{\beta_1} \dots \delta^{\mu_n}_{\alpha_n} \delta^{\nu_n}_{\beta_n]} \prod_{k=1}^{n} R^{\alpha_k \beta_k}_{\mu_k \nu_k}. \qquad (8.1)
$$

The definition of the antisymmetric product and relations to analogue expressions in terms of
Levi-Civita symbols is given in appendix A. In spite of being higher-derivative operators, the
LL invariants are believed to avoid inconsistencies that are usually introduced by such terms
This is due to the antisymmetric structure in (8.1) which ensures the absence of more than two
time derivatives acting on one field in the action. It is interesting that the same structure
appears in the interaction potential of ghost-free bimetric theory, but of course this is not a
complete coincidence: in the decoupling limit of nonlinear massive gravity it can be seen that
the antisymmetric structure removes higher time derivatives from the longitudinal modes of
the Stückelberg fields [57, 58].

8.1.1. *Extended action and equations of motion.* In order to extend ghost-free bimetric theory to d-dimensions, we first write the Hassan–Rosen bimetric action in the more general form

$$
\begin{equation}
\begin{split}
S_{\text{HR}} ={}& m_g^{d-2} \int d^d x \sqrt{g} R(g) + m_f^{d-2} \int d^d x \sqrt{f} R(f) \\
& - 2m^d \int d^d x \sqrt{g} \sum_{n=0}^{d} \beta_n e_n (\sqrt{g^{-1}f}),
\end{split}
\tag{8.2}
\end{equation}
$$

where the elementary symmetric polynomials are still defined through the same recursion
formula (A.2). The consistency proof of [76], as outlined in sections 3.4 and 5.2, generalises
straightforwardly to the higher-dimensional case. It has also been shown that the mass
spectrum is very similar to the four-dimensional case and contains a massless and massive
excitation [143].

Next, guided by the above considerations, we extend the Hassan–Rosen action by the
following terms, as done in [143, 268]

$$
S_{\text{LL}} = m_g^{d-2} \int d^d x \sum_{n=2}^{[d/2]} \sqrt{-g} l_n^g \mathcal{L}_{(n)}(g) + m_f^{d-2} \int d^d x \sum_{n=2}^{[d/2]} \sqrt{-f} l_n^f \mathcal{L}_{(n)}(f), \quad (8.3)
$$

where we have introduced two sets of couplings, $l_n^g$ and $l_n^f$, of mass dimension 2(1 − n). It is
well known that the LL invariant $\mathcal{L}_{(n)}$ is topological in $d = 2n$ and vanishes for $2n > d$ which
is why the sums in (8.3) terminate at integer part of $d/2$.

The equations of motion for g<sub>μν</sub> and f<sub>μν</sub> that follow from the LL extended bimetric action
are

<footnote n="44">⁴⁴ See, however, [260–262]. Moreover, note that the structure of the kinetic terms is of course only determined up to field redefinitions.</footnote>
$$G_{\mu\nu}(g) + \sum_{n=2}^{[d/2]} l_n^g G_{\mu\nu}^{(n)}(g) + \frac{m^d}{m_g^{d-2}} V_{\mu\nu}^g = 0,$$

$$G_{\mu\nu}(f) + \sum_{n=2}^{[d/2]} l_n^f G_{\mu\nu}^{(n)}(f) + \frac{m^d}{m_g^{d-2}} V_{\mu\nu}^f = 0. \qquad (8.4)$$

Here, $G_{\mu\nu}^{(n)}$ are the Lovelock tensors that follow from variation of (8.3) and $V_{\mu\nu}^g$ and $V_{\mu\nu}^f$ are the contributions from the interaction potential which were defined in (5.3) and in which the summations now run from 0 to d. We will not need the explicit expressions for the Lovelock tensors, but let us note the important property

$$g^{\mu\nu} G_{\mu\nu}^{(n)}(g) = \frac{2n-d}{2} \mathcal{L}_{(n)}(g), \qquad (8.5)$$

which will be useful later on.

### 8.1.2. Proportional backgrounds.
As in pure bimetric theory in $d=4$ we are interested in finding the proportional background solutions to (8.4) and study perturbations around those. Again, these backgrounds correspond to maximally symmetric spacetimes, for which the curvatures satisfy

$$R_{\mu\nu\rho\sigma}(\bar{g}) = \frac{2\lambda}{(d-1)(d-2)}(\bar{g}_{\mu\rho}\bar{g}_{\nu\sigma} - \bar{g}_{\nu\rho}\bar{g}_{\mu\sigma}), \\ R_{\mu\nu}(\bar{g}) = \frac{2\lambda}{d-2}\bar{g}_{\mu\nu}, \qquad R(\bar{g}) = \frac{2d\lambda}{d-2}, \qquad (8.6)$$

with cosmological constant $\lambda$. On such backgrounds with constant curvature, the nth order LL invariant is proportional to $\lambda^n$

$$\mathcal{L}_{(n)} = N_n(d) \lambda^n, \quad \text{with} \quad N_n(d) = \frac{2^n d!}{(d-1)^n (d-2)^n (d-2n)!} \qquad (8.7)$$

Next, we trace the equations of motion (8.4) in order to obtain two scalar equations and insert (8.6) for the curvatures. Furthermore, we make use of the identity (8.5) and plug in the proportional ansatz, $\bar{f}_{\mu\nu} = c^2 \bar{g}_{\mu\nu}$. Finally, in the equation for $f_{\mu\nu}$, we use $\mathcal{L}_{(n)}(c^2 \bar{g}) = c^{-2n} \mathcal{L}_{(n)}(\bar{g})$. The resulting two equations read

$$\lambda + \sum_{n=2}^{[d/2]} l_n^g \frac{d-2n}{2d} N_n(d) \lambda^n - \Lambda_g = 0, \\ \lambda + \sum_{n=2}^{[d/2]} c^{2-2n} l_n^f \frac{d-2n}{2d} N_n(d) \lambda^n - \Lambda_f = 0. \qquad (8.8)$$

Here, $\Lambda_g$ and $\Lambda_f$ are the contributions from the potential evaluated on the proportional backgrounds which were defined earlier in (5.11). Note, however, that these no longer correspond to the physical cosmological constants and are not necessarily equal.

For general parameters $\beta_n$, $l_n^g$, and $l_n^f$, the two algebraic equations in (8.8) serve to determine the background curvature $\lambda$ as well as the proportionality constant $c^2$. Note also that in the absence of the LL contributions, i.e. for $l_n^g = l_n^f = 0$, we re-arrive at the bimetric background equations, $R(\bar{g}) = \frac{2d}{d-2} \Lambda_g$ and $\Lambda_g = \Lambda_f$.

### 8.1.3. Mass spectrum.
Computing the linearised equations of motion around general background solutions to (8.4) is quite involved, even when the potential contributions vanish.
However, as has been observed in [267], around constant curvature backgrounds with cosmological constant $\lambda$, the fluctuation equations of GR augmented by the LL terms with couplings $l_n^g$ assume the same form as the ones obtained from Einstein's equations. The only difference is that, in the linearised equations of the LL extended theory, the Planck mass $m_g$ needs to be replaced by an effective mass parameter $\tilde{m}_g$ defined through

$$ \frac{\tilde{m}_g^{d-2}}{m_g^{d-2}} = 1 + (d-3)! \sum_{n=2}^{[d/2]} \frac{n(d-2n)l_n^g}{(d-2n)!} \left( \frac{2\lambda}{(d-1)(d-2)} \right)^{n-1} . \quad (8.9) $$

It is a remarkable feature of the Lovelock invariants that, in spite of being higher-derivative operators in the nonlinear theory, they give rise to linearised equations that are only second order in derivatives.

For the bimetric case the above result implies that the linear equations for the theory with LL terms are the same as in pure bimetric theory if one replaces $m_g$ by $\tilde{m}_g$ as above and $m_f$ by $\tilde{m}_f$ given by

$$ \frac{\tilde{m}_f^{d-2}}{m_f^{d-2}} = 1 + (d-3)! \sum_{n=2}^{[d/2]} \frac{n(d-2n)c^{2-2n}l_n^f}{(d-2n)!} \left( \frac{2\lambda}{(d-1)(d-2)} \right)^{n-1} . \quad (8.10) $$

Hence, using the results of [267], without any further computation we can conclude that the linear equations in a maximally symmetric background with cosmological constant $\lambda$ will again diagonalise into a massless and massive equation, where in the latter the Fierz-Pauli mass (5.21) is now given in terms of $\tilde{m}_g$ and $\tilde{m}_f$.

### 8.1.4. Partial masslessness in $d > 4$
To conclude the discussion of the LL extension of bimetric theory, let us comment on the possibility to extend the linear symmetry of PM fields to the nonlinear level within this generalised setup. Since the proportional backgrounds and their perturbations possess the same structure as in pure bimetric theory in $d = 4$, the linear equations are again invariant under a PM symmetry if the mass is on the Higuchi bound which in $d$-dimensions corresponds to the value $m_{FP}^2 = \frac{2}{d-1}\lambda$. The arguments given in section 7.2.1 leading to a PM candidate theory can straightforwardly be applied to the higher-dimensional case. The criterion of leaving the proportionality constant $c$ undetermined by the background equations (8.8) singles out unique PM candidate theories which have been identified in [143] up to $d = 8$. All the interaction parameters are fixed with respect to each other in these models and, in particular, the LL coefficients $l_n^g$ and $l_n^f$ are nonzero. This result confirms earlier findings [228] which revealed that the linear PM symmetry cannot be extended in $d > 4$ unless higher-derivative terms are included into the action.

## 8.2. Multiple interacting spin-2 fields

Ghost-free bimetric theory contains the correct number of degrees of freedom for a massless and a massive spin-2 mode. While no-go theorems forbid consistent interactions among more than one massless spin-2 field [258], it is possible to extend bimetric theory by additional massive spin-2 degrees of freedom. Hinterbichler and Rosen first wrote down these ghost-free interactions using vierbeins instead of metrics<fnref n="45" /> [269]. In the following we will derive the consistent multi-spin-2 interactions in the metric formulation and then review their formulation in terms of vierbeins. More work on the vierbein formulation of massive gravity, bimetric and multimetric theory can be found in [68, 101, 102, 136, 270–282].

<footnote n="45">⁴⁵ The formulation of the full ghost-free bimetric potential in terms of vierbeins first appeared in [153].</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1. The centre couplings (left) and the chain couplings (right).</figcaption></figure>

8.2.1. Multiple bimetric couplings. When constructing a theory for more than two spin-2 fields, of course, the main requirement on its structure remains the absence of the Boulware–Deser ghosts. The interaction potential for $\mathcal{N}$ fields therefore has to be chosen such that the action in Hamiltonian formulation gives rise to $(4 + \mathcal{N} - 1)$ constraints. Out of these, four should be first class constraints, generating the group of GCT, while the remaining $(\mathcal{N} - 1)$ are expected to be second class and serve to remove all Boulware–Deser ghosts from the spectrum.

An obvious way to satisfy these requirements is to simply add up several copies of the Hassan–Rosen action (5.1) for fields $f_{\mu\nu}$ and $g_{\mu\nu}(I)$ with $I = 1, \dots, \mathcal{N} - 1$, to obtain a theory in which one of the metrics is coupled to all the others

$$S_{\text{centre}}[f, \{g(I)\}] = \sum_{I=1}^{\mathcal{N}-1} S_{\text{HR}}[g(I), f]. \qquad (8.11)$$

This will be consistent because in every interaction potential one can redefine the ADM variables for $g_{\mu\nu}(I)$ to make it linear in the lapse and shift of $f_{\mu\nu}$ as well as the lapse of $g_{\mu\nu}(I)$, such that one obtains the desired $(4 + \mathcal{N} - 1)$ constraints.

The left panel of figure 1 shows how such a coupling can be visualised in a graph [269]. Each solid dot represents a different spin-2 field and the lines stand for standard Hassan–Rosen bimetric couplings between them. Fields corresponding to dots that are not connected by a line do not directly interact. Due to the corresponding picture with one metric in the centre, we shall refer to interactions of the form (8.11) as ‘centre couplings’.

Another option to obtain consistent multiple bimetric interactions is to build a ‘chain’ of $\mathcal{N}$ coupled spin-2 fields $g_{\mu\nu}(I)$, where $I = 1, \dots, \mathcal{N}$. That is to say, we take an action of the form

$$S_{\text{chain}}[\{g(I)\}] = \sum_{I=1}^{\mathcal{N}-1} S_{\text{HR}}[g(I), g(I+1)], \qquad (8.12)$$

for which the corresponding graph is depicted in the right panel of figure 1. This construction works because in the first term one can redefine the ADM variables for $g(1)$ such that it becomes linear in the lapse and shift of $g(2)$ as well as the lapse of $g(1)$. In the second term one performs a similar redefinition for the shift of $g(2)$ which, since the redefinition is linear in the lapses, does not destroy the linearity in the lapse of $g(2)$ in the first term. Note however that the second term will no longer be linear in the shift of $g(2)$. Continuing with this procedure in all the terms results in an action that is linear in the lapse and shift of $g(\mathcal{N})$ and in all the other lapses.

In fact, it is also possible to combine the two above constructions in the following manner: first introduce a chain coupling for fields $f(K)$, $K = 1, \dots, K$, and then to each $f(K)$ attach $\mathcal{I}_K$ fields $g_{(K)}(I_K)$, $I_K = 1, \dots, \mathcal{I}_K$, in a centre coupling. This gives rise to the most
<figure><img src="image_6.png" /><figcaption>**Figure 2.** The combined chain and centre couplings (**left**) and a forbidden loop coupling (**right**).</figcaption></figure>

general ghost-free multimetric action

$$S_{\text{multi}} = S_{\text{chain}}[\{f(K)\}] + \sum_{K=1}^{\mathcal{K}} S_{\text{centre}}[f(K), \{g_{(K)}(I_K)\}]. \quad (8.13)$$

In this theory the total number of fields is $\mathcal{N} = \mathcal{K} + \sum_{K=1}^{\mathcal{K}} \mathcal{I}_k$. An example for a graph is shown in the left panel of figure 2. Note that the action $S_{\text{multi}}$ still consists purely of multiple copies of the Hassan–Rosen bimetric action; we have not introduced any type of new interaction term. Of course, in the above construction, it is also possible to add further chain or centre couplings to each of the $g_{(K)}(I_K)$. However, it is important to realise that not all possible bimetric couplings are allowed. In particular, one may not couple the metrics in a ‘loop’, which means adding terms such as $S_{\text{HR}}[f(K), f(K+2)]$ or $S_{\text{HR}}[f(K+1), g_{(K)}(I_K)]$ to (8.13). To see why this fails, consider a simple loop coupling among three metrics

$$S_{\text{loop}}(g, f, h) = S_{\text{HR}}(g, f) + S_{\text{HR}}(f, h) + S_{\text{HR}}(h, g), \quad (8.14)$$

whose graph is depicted in the right panel of figure 2. Now imagine that in the first term one redefines the shift of $g$ to render it linear in the lapse and shift of $f$ and the lapse of $g$. Then, in the second term redefine the shift of $f$ such that the first two terms afterwards are linear in the lapse and shift of $h$ and the lapses of $g$ and $f$. In the third term, one would now have to redefine the shift of $h$, but from the particular form (3.19) of the redefinition it is clear that afterwards, the action will again contain the original shift variable for $g$ which cannot be expressed in terms of the redefined variables only. Therefore there is at least no obvious consistent redefinition of ADM variables that renders (8.14) linear in all the lapses and we conclude that the loop couplings most likely give rise to Boulware–Deser ghosts. This is in agreement with results of a similar analysis performed in [283].

### 8.2.2. Vierbein formulation.
The multimetric action (8.13) consists of several copies of the Hassan–Rosen action (5.1). It is possible to rewrite the latter in terms of different variables which avoid the appearance of the square-root matrix$^{46}$. To this end, let us introduce the tetrads or vierbeins defined for each metric

$$g_{\mu\nu}(I) = E_{\mu}^{a}(I) \eta_{ab} E_{\nu}^{b}(I), \quad (8.15)$$

or in matrix notation, $g(I) = E^T(I)\eta E(I)$. The expression for the metric in terms of vierbeins is invariant under local Lorentz transformations $\Lambda_b^a(I)$ that act on the upper index of $E_{\mu}^{b}(I)$ according to $E(I) \to \Lambda(I)E(I)$ with $\Lambda^T\eta\Lambda = \eta$. Associated to each vierbein, there are thus six gauge invariances which serve to reduce the number of physically relevant components from 16 to 10. Any metric theory that is reformulated in terms of vierbeins is automatically

<footnote>$^{46}$ For suggestions to remove the square-root matrix by introducing auxiliary fields, see [64, 284].</footnote>
invariant under these transformations and, in particular, the Einstein–Hilbert terms $S_{EH}(I)$, in which the metric $g(I)$ is replaced by the vierbein, are Lorentz invariant.

In terms of the respective vierbeins, the square-root matrix $\sqrt{g(1)^{-1}g(2)}$ built out of two metrics $g(1)$ and $g(2)$ becomes

$$S = \sqrt{g(1)^{-1}g(2)} = \sqrt{E(1)^{-1}\eta^{-1}(E(1)^{\mathrm{T})^{-1}E(2)^{\mathrm{T}}\eta E(2)}}, \quad (8.16)$$

It is now convenient to fix the Lorentz gauge for one of the vierbeins such that $E(2)E(1)^{-1}\eta^{-1}$ is a symmetric matrix. We thus impose the so-called ‘symmetric’ gauge (or Deser-van Nieuwenhuizen gauge [285])

$$\eta^{-1}(E(1)^{\mathrm{T}})^{-1}E(2)^{\mathrm{T}} = E(2)E(1)^{-1}\eta^{-1}. \quad (8.17)$$

In this particular Lorentz frame, we can easily evaluate the square-root as

$$S = \sqrt{E(1)^{-1}E(2)E(1)^{-1}E(2)} = E(1)^{-1}E(2). \quad (8.18)$$

Using the identity

$$\epsilon_{c_1 \dots c_d} = (\det E) \epsilon_{\nu_1 \nu_2 \nu_3 \nu_4} (E^{-1})^{\nu_1}_{\quad c_1} (E^{-1})^{\nu_2}_{\quad c_2} (E^{-1})^{\nu_3}_{\quad c_3} (E^{-1})^{\nu_4}_{\quad c_4}, \quad (8.19)$$

we can then write the bimetric action in the form

$$S_{\mathrm{HR}} = S_{\mathrm{EH}}(1) + S_{\mathrm{EH}}(2) - S_{\mathrm{int}}, \quad S_{\mathrm{int}} \equiv 2m^4 \int d^4x \sum_{n=0}^{4} V_n, \quad (8.20)$$

where $S_{EH}(I)$ are the standard Einstein–Hilbert actions for $E(I)$ and

$$V_n = b_n \epsilon^{\mu_1 \dots \mu_n \mu_{n+1} \dots \mu_4} \epsilon_{a_1 \dots a_n a_{n+1} \dots a_4} E^{a_1}_{\quad \mu_1} (1) \dots E^{a_n}_{\quad \mu_n} (1) E^{a_{n+1}}_{\quad \mu_{n+1}} (2) E^{a_4}_{\quad \mu_4} (2), \quad (8.21)$$

where $b_n = \frac{2\beta_n}{n!(4-n)!}$ are the same parameters as in (3.33). Since we have only fixed one of the two Lorentz gauges, this interaction potential still has one overall local Lorentz invariance under which both vierbeins transform in the same way.

It is important to note that the vierbeins in (8.21) must satisfy the symmetrisation condition (8.17), otherwise the vierbein formulation is not equivalent to bimetric theory for general values of $b_n$. More precisely, if we started from (8.21) without imposing the condition by hand, then the dynamics of the vierbein theory would allow for configurations that do not give back the ghost-free bimetric formulation. As discussed in [269], for each vielbein there exist six combinations of its equations of motion that do not contain any derivatives. These read

$$\frac{\delta S_{\mathrm{int}}}{\delta E^{a}_{\mu}(I)} E^{c}_{\mu}(I) \eta_{cb} - \frac{\delta S_{\mathrm{int}}}{\delta E^{b}_{\mu}(I)} E^{c}_{\mu}(I) \eta_{ca} = 0. \quad (8.22)$$

For models with $b_2 = b_3 = 0$ or $b_1 = b_2 = 0$, these equations imply the symmetrisation condition (8.17) and the equivalence to bimetric theory is dynamically guaranteed. For more general parameters, the symmetrised vierbeins still solve the equations but other solutions also exist [286, 287]. These disconnected branches give rise to different theories which contain the Boulware–Deser ghost instability [288].

Another important aspect, first pointed out in [286], is that one cannot always arrive at the symmetrised form (8.17) by acting with a Lorentz transformation on general vierbeins. For configurations that do not allow this, the square-root cannot be evaluated as in (8.18). In fact, the condition on the combination of vierbeins to be symmetrisable by a Lorentz transformation is exactly the same as the one we derived in section 3.3.3 on the metric components by requiring that the square-root exists [99]: ensuring that the Lorentz boost velocity $v_a$ satisfies $v_a \delta^{ab}v_b < 1$ is equivalent to requiring the existence of a real solution to the
<figure><img src="image_7.png" /><figcaption>Figure 3. A three- and a four-point vertex.</figcaption></figure>

scalar square-root $\sqrt{x}$ that appears in the ADM decomposition of $\sqrt{g^{-1}f}$ in (3.27). This condition in turn ensures that one can always impose the ‘symmetric’ gauge by a Lorentz transformation.

### 8.2.3. More general interactions?
So far we have used the earlier results of the bimetric ADM analysis in order to couple several spin-2 fields to each other using only bimetric interactions. The next question is whether there exist more general couplings than the bimetric vertex. We introduce graphs for two examples of such n-point vertices in figure 3. A priori, such couplings could be very different from the bimetric ones and the ADM analysis has to be redone from the start.

In [269] it was proposed that the interactions of vierbeins could be generalised to

$$ S_{\text{int}} = \frac{m^4}{4} \sum_{\{I\}} \int d^4 x \, T^{I_1 I_2 I_3 I_4} U_{I_1 I_2 I_3 I_4}, $$ 
with $U_{I_1...I_4} = \epsilon^{\mu_1 \mu_2 \mu_3 \mu_4} \epsilon_{a_1 a_2 a_3 a_4} E^{a_1}_{\mu_1}(I_1) E^{a_2}_{\mu_2}(I_2) E^{a_3}_{\mu_3}(I_3) E^{a_4}_{\mu_4}(I_4)$, (8.23)

where $T^{I_1 I_2 I_3 I_4}$ are coupling constants symmetric in all the indices $I_n \in \{I\}$ which run from 1 to $N$. The total multivielbein action would then be given by

$$ S_{\text{multi}} = \sum_{I=1}^{N} S_{\text{EH}}(I) + S_{\text{int}}. \qquad (8.24) $$

Note that the antisymmetric structure of (8.23) ensures that there cannot be more than four vierbeins (or in d-dimensions, not more than d vielbeins) interacting in one vertex.

However, the results obtained in [288] severely constrain the structure of the ‘tensor’ of coupling constants $T^{I_1...I_4}$ and only diagrams of the type displayed in figure 1 are allowed. The more general couplings displayed in figure 3 which would arise from an unconstrained $T^{I_1...I_4}$ and give rise to diagrams with at most four lines ending in one vertex are excluded. The same holds for the loop couplings in the right panel of figure 2. Although these forbidden terms have been shown to possess a ghost-free decoupling limit, they do spoil the consistency of the full theory [136, 288] and the reason can be traced back to the fact that the corresponding equations of motion are incompatible with the symmetrisation condition (8.17). In summary, the only consistent vierbein couplings are those that give rise to a metric formulation of the type (8.13).

## 9. Discussion and outlook

In this article we summarised the recent developments of nonlinear theories involving massive spin-2 fields. It was shown how to construct the unique ghost-free interaction potential of nonlinear massive gravity with fixed reference metric, which we afterwards generalised to the
fully dynamical bimetric theory. The latter describes the interactions of a massless spin-2 field with a massive one and, when coupled to matter, captures the behaviour of gravity in the presence of an additional tensor field. We have seen that there exists a limit in the parameter space of bimetric theory which takes the theory arbitrarily close to GR. In this limit, viable cosmological solutions can be found where, in the absence of vacuum energy, the dark energy scale is set by the spin-2 mass which is protected against receiving large quantum corrections. We have discussed the idea of nonlinear theories for PM spin-2 fields and in the process established a connection between bimetric theory and conformal gravity. Finally, we have reviewed the extension of bimetric theory to dimensions greater than four and to couplings between more than two spin-2 fields.

On the phenomenological side, one of the most interesting open questions is what possible implications the existence of a massive spin-2 field in nature could have on the dark matter problem. Several scenarios where dark matter components are attributed to the matter sector of the second metric $f_{\mu\nu}$ have already been considered in [149–152]. From our point of view, however, the more interesting option is to regard the second metric as the dark matter field. Since the spin of the dark matter particle is unknown, such a scenario is not immediately excluded and, as we already discussed in section 5.4.2, the largeness of the physical Planck mass automatically weakens the interactions of the massive spin-2 field with all Standard Model particles.

On the theoretical side, it would be of great importance to understand the origin of the mass term for the spin-2 field at a fundamental level. For spin-1 fields, it is well-known that the consistent picture at the quantum level is to generate their mass via spontaneous breaking of the gauge symmetry in the massless theory. To this end, it is necessary to invoke additional (Higgs) fields whose vacuum expectation values break gauge invariance and set the spin-1 mass scale. Since aspects of field theories generically do not become simpler with increasing spin, it is natural to expect an underlying mechanism that generates the mass term also for spin-2 fields. Such a ‘spin-2 Higgs mechanism’ should break the two independent copies of diffeomorphism invariance in the $g_{\mu\nu}$ and $f_{\mu\nu}$ sectors down to their diagonal subgroup, which is the residual gauge symmetry of bimetric theory. The precise nature of the fields that are needed to trigger this symmetry breaking spontaneously is still unknown; a realisation on anti-de-Sitter space has been suggested in [289, 290]. For a more detailed discussion of the topic, we refer the reader to the recent article [291].

Related to the origin of the mass term is the behaviour of bimetric theory at the quantum level. Quantum corrections at one-loop order in the effective field theory picture have been computed in massive gravity [214] and bimetric theory [170] and the results confirmed the stability of the spin-2 mass scale. On the other hand, without knowing the underlying mechanism that generates the interaction potential for the two metrics, it is difficult to tell to what extend these results can be trusted.

Theories with extended symmetries are known to generically exhibit an improved quantum behaviour. An example for the spin-2 case is the action for conformal gravity with its Weyl symmetry. Unfortunately, this theory suffers from a ghost but we have seen that bimetric theory contains a ghost-free model which seems to be closely related to the Weyl invariant action. Understanding its properties further and investigating possibilities to realise the gauge symmetry for PM at the nonlinear level (most likely by invoking additional degrees of freedom) could give us important new insights on the nature of quantum gravity.

Since consistent theories involving spin-2 fields are so rare, we do not expect to encounter a large variety of possibilities to extend them further. As it stands now, it does not seem to be possible to extend the kinetic sector in four-dimensions [135, 136, 259] and the interaction potential is the only consistent one by construction. Up to field redefinitions, the
structure of the ghost-free bimetric action is therefore unique. On the other hand, it is of course worth studying couplings to fields with different and, in particular, higher spins. The spectrum of presently known higher-spin theories contains only massless spin-2 fields (see, for instance, [292–294]). It would be interesting to see whether bimetric theory can give hints on the form of more general interactions including massive spin-2 degrees of freedom. At this stage, we can only hope that one day we will be able to write down the most general interactions for massless and massive spin-2 fields.

## Acknowledgments

We would like to thank Yashar Akrami, Eugeny Babichev, Laura Bernard, Cédric Deffayet, Jonas Enander, Dario Francia, Aleksander Garus, Fawad Hassan, Frank Könnig, Alessandro Sfondrini and Adam Solomon for fruitful discussions and comments on the draft. We are particularly grateful to Niklas Beisert who motivated us to write this review.

The research of MvS leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/ 2007-2013 Grant Agreement no. 307934). ASM is supported by ERC grant no. 615203 under the FP7 and the Swiss National Science Foundation through the NCCR SwissMAP.

## Appendix A. Technical details of interaction potential

In the text, for notational simplicity, we frequently write the interactions in terms of elementary symmetric polynomials $e_n(S)$. These can be defined in various equivalent ways, for example via the following commonly used equalities which hold for any square $d \times d$ matrix $S$,

$$ 
\begin{aligned}
e_n(S) &= S_{[\mu_1}^{\mu_1} \cdots S_{\mu_n]}^{\mu_n} = \frac{1}{n!} \delta_{\nu_1 \cdots \nu_n}^{\mu_1 \cdots \mu_n} S_{\mu_1}^{\nu_1} \cdots S_{\mu_n}^{\nu_n} \\
&= \frac{1}{n!(d-n)!} \epsilon^{\mu_1\mu_2\cdots\mu_n\lambda_{n+1}\cdots\lambda_d} \epsilon_{\nu_1\nu_2\cdots\nu_n\lambda_{n+1}\cdots\lambda_d} S_{\mu_1}^{\nu_1} \cdots S_{\mu_n}^{\nu_n}.
\end{aligned} 
\quad (A.1) 
$$

Here we have expressed them in terms of antisymmetrisation, the generalised Kronecker delta and the Levi-Civita symbol respectively, all normalised by unit weight. Obviously the last equality only makes sense when $d$ is the dimension spanned by the indices. They can also be defined through a dimension-independent recursion formula

$$ e_n(S) = \frac{(-1)^{n+1}}{n} \sum_{k=0}^{n-1} (-1)^k \mathrm{Tr}(S^{n-k}) e_k(S), \quad \text{with } e_0(S) = 1. \quad (A.2) $$

In $d$-dimensions, they satisfy $e_d(S) = \det(S)$ and $e_n(S) = 0$ for $n > d$. Their explicit expressions for $n \le 4$ read

$$ 
\begin{aligned}
e_0(S) &= 1, & e_1(S) &= [S], & e_2(S) &= \frac{1}{2}([S]^2 - [S^2]), \\
e_3(S) &= \frac{1}{6}([S]^3 - 3[S][S^2] + 2[S^3]) & & \\
e_4(S) &= \frac{1}{24}([S]^4 - 6[S]^2[S^2] + 3[S^2]^2 + 8[S][S^3] - 6[S^4]), & &
\end{aligned} 
\quad (A.3) 
$$

where square brackets denote a matrix trace. Moreover, the $e_n$ obey the following useful identity for any matrix $S$ and parameter $\lambda$
$$e_n(1 + \lambda S) = \sum_{k=0}^{n} \binom{d-k}{n-k} \lambda^k e_k(S). \qquad (A.4)$$

A special case ($n = d$) of this is

$$\det(1 + \lambda S) = \sum_{n=0}^{d} \lambda^n e_n(S). \qquad (A.5)$$

This means that the consistent interaction potential can be regarded as a *deformed determinant* [46]. The variation of the $e_n(S)$ can be computed from e.g. (A.2) and is given by$^{47}$

$$\delta e_n(S) = - \sum_{k=1}^{n} (-1)^k \mathrm{Tr}[S^{k-1} \delta S] e_{n-k}(S), \quad n \ge 1, \qquad (A.6)$$

with $\delta e_0(S) = 0$ (since $e_0(S) = 1$). Together with $2\mathrm{Tr}[S^{k-1}\delta S] = \mathrm{Tr}[S^{k-2}\delta S^2]$ and $\delta S^2 = -g^{-1}\delta g S^2 + g^{-1}\delta f$ (which follows from $S^2 = g^{-1}f$) it is then straightforward to derive the field equations of the theory. Explicitly, the contributions from the interaction potential to the equations of motion in (5.3) are

$$(Y_{(n)})^{\mu}_{\nu}(S) \equiv \sum_{k=0}^{n} (-1)^k (S^{n-k})^{\mu}_{\nu} e_k(S). \qquad (A.7)$$

In $d$-dimensions, they satisfy $Y_{(n)} = 0$ for $n \ge d$. In particular, $Y_{(d)} = 0$ is simply the statement of the Cayley–Hamilton theorem, that any square matrix satisfies its own characteristic (or secular) equation.

Finally we note that all of the above expressions can easily be rewritten in terms of the inverse $S^{-1}$ by using the identity

$$e_n(S^{-1}) = \frac{e_{d-n}(S)}{e_d(S)}. \qquad (A.8)$$

## Appendix B. Derivation of redefined shift vector

Here we derive the expressions for the redefinition (3.19) of the ADM shift vector $N^i = \gamma^{ij}N_j$ and the matrices $A$ and $B$ in (3.27), following [61]. We start from the most general ansatz for the redefinition

$$N^i = c^i + Nd^i, \qquad (B.1)$$

where $c^i$ and $d^i$ are functions of $\gamma_{ij}$, the new shift vector $n^i$ and the ADM components of $f_{\mu\nu}$. They will be determined in what follows. Recall that the redefinition has to be linear in the lapse $N$ in order not to introduce nonlinearities into the Einstein–Hilbert term (3.4). Next we turn to equation (3.18) and take the square of both sides to arrive at

$$g^{-1}f = \frac{1}{N^2}A^2 + \frac{1}{N}(AB + BA) + B^2. \qquad (B.2)$$

The left-hand side of this equation can be expressed in terms of ADM variables using (3.3) and (3.12). We get

<footnote>$^{47}$ See e.g. [105] for an explicit derivation.</footnote>
$$g^{\mu\rho} f_{\rho\nu} = \frac{1}{N^2} \left( \begin{array}{cc} L^2 - L_l \phi^{lk} L_k + N_l \gamma^{lk} L_k & -L_j + N_l \gamma^{lk} \phi_{kj} \\ N^2 \gamma^{ik} L_k - N^i (L^2 - L_l \phi^{lk} L_k + N_l \gamma^{lk} L_k) & N^2 \gamma^{ik} \phi_{kj} + N^i (N^k \phi_{kj} - L_j) \end{array} \right). \qquad (B.3)$$

In this expression we replace $N^i$ in terms of $n^i$ using our ansatz (B.1). Next, we collect the pieces according to their order in $1/N$, which gives

$$g^{-1}f = \frac{1}{N^2} \mathbb{E}_0 + \frac{1}{N} \mathbb{E}_1 + \mathbb{E}_2, \qquad (B.4)$$

where, in terms of $a_0 \equiv L^2 - L_l \phi^{lk} L_k + c^l L_l$ and $a_i \equiv -L_i + c^i \phi_{li}$, we have defined the matrices

$$\mathbb{E}_0 = \begin{pmatrix} a_0 & a_j \\ -a_0 c^i & -c^i a_j \end{pmatrix}, \quad \mathbb{E}_1 = \begin{pmatrix} d^l L_l & d^l \phi_{lj} \\ -(d^l L_l c^i + a_0 d^i) & -(c^i d^l \phi_{lj} + d^i a_j) \end{pmatrix} \\ \mathbb{E}_2 = \begin{pmatrix} 0 & 0 \\ (\gamma^{il} - d^i d^l)L_l & (\gamma^{il} - d^i d^l)\phi_{lj} \end{pmatrix}. \qquad (B.5)$$

Comparing this to the right-hand side of (B.2) and equating the coefficients in front of equal powers of $1/N$, we see that

$$\mathbb{A}^2 = \mathbb{E}_0, \quad \mathbb{B}^2 = \mathbb{E}_2, \quad \mathbb{A}\mathbb{B} + \mathbb{B}\mathbb{A} = \mathbb{E}_1. \qquad (B.6)$$

Since $\mathbb{E}_0$ is a projector on a one-dimensional subspace, its square-root can easily be evaluated, resulting in

$$\mathbb{A} = \sqrt{\mathbb{E}_0} = \frac{1}{L\sqrt{x}} \begin{pmatrix} a_0 & a_j \\ -a_0 c^i & -c^i a_j \end{pmatrix}, \quad x \equiv \frac{1}{L^2} (a_0 - c^i a_j). \qquad (B.7)$$

The special structure of $\mathbb{E}_2$ also allows us to evaluate its square-root which gives

$$\mathbb{B} = \sqrt{\mathbb{E}_2} = \sqrt{x} \begin{pmatrix} 0 & 0 \\ D_k^i L^k & D_j^i \end{pmatrix}, \quad D_j^i \equiv \frac{1}{\sqrt{x}} \sqrt{(\gamma^{il} - d^i d^l)\phi_{lj}}. \qquad (B.8)$$

Using these expressions for $\mathbb{A}$ and $\mathbb{B}$ in the third equation of (B.6), we obtain

$$\mathbb{A}\mathbb{B} + \mathbb{B}\mathbb{A} = \mathbb{E}_1 \quad \Rightarrow \quad d^i = \frac{1}{L} D_k^i (c^k - \phi^{kl} L_l), \qquad (B.9)$$

where we have also made use of the symmetry property $\phi_{ik} D_j^k = \phi_{jk} D_k^i$ which follows directly from the definition of $D$ in (B.8). This equation determines a particular combination of the unknown functions $d^i$ and $c^i$ in the redefinition. The redefinition is thus not uniquely determined. One simple option to fix the ambiguity is to choose the new shift vectors as $n^k = \frac{1}{L}(c^k - \phi^{kl} L_l)$. Then the above condition reduces to $d^i = D_k^i n^k$ and the redefinition becomes

$$N^i = c^i + Nd^i = L^i + Ln^i + ND_k^i n^k, \qquad (B.10)$$

where $L^i \equiv \phi^{ij} L_j$. The definition of the matrix $D$ in (B.8) depends on $d^i$ and therefore gives rise to the following matrix equation that needs to be solved for $D$

$$\sqrt{x} D = \sqrt{(\gamma^{-1} - Dn(Dn)^T)\phi}, \quad x \equiv 1 - n^l \phi_{lk} n^k. \qquad (B.11)$$

As shown in [61], the solution to this equation is given by (3.20).
## Appendix C. Short summary of standard GR cosmology

Here we very briefly recapitulate the derivation of the cosmological evolution equations in GR. Consider Einstein's equations of motions for $g_{\mu\nu}$

$$R_{\mu\nu} - \frac{1}{2}g_{\mu\nu}R = \frac{1}{M_{\text{Pl}}^2}(T_{\mu\nu} - \rho_{\Lambda}g_{\mu\nu}), \qquad (C.1)$$

where $T_{\mu\nu}$ is the stress-energy tensor obtained from variation of the matter Lagrangian and we have defined the constant energy density $\rho_{\Lambda} \equiv \Lambda M_{\text{Pl}}^2$. When looking for homogeneous and isotropic solutions in GR, one makes a Friedmann–Robertson–Walker ansatz for the metric

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

in which $k = 0, -1, +1$ parameterises the curvature of the Universe (flat, open, closed), $a(t)$ is the scale factor, and $d\Omega = d\theta^2 + \sin^2\theta d\phi^2$. In accordance with homogeneity and isotropy, the stress-energy tensor is taken to be that of a perfect fluid

$$(T^g)^{\mu}_{\nu} = \text{diag}(-\rho, p, p, p), \qquad (C.3)$$

where $\rho(t)$ and $p(t)$ are the time dependent energy density and pressure of the fluid, respectively. As a consequence of the Bianchi identity, the source is automatically covariantly conserved which implies the continuity equation

$$\dot{\rho} + 3\frac{\dot{a}}{a}(\rho + p) = 0, \qquad (C.4)$$

where dots denote time derivatives. Plugging the diagonal ansätze (C.2) and (C.3) into Einstein's equations (C.1) gives a set of four equations. The 00-component is Friedmann's equation

$$\left(\frac{\dot{a}}{a}\right)^2 + \frac{k}{a^2} = \frac{\rho + \rho_{\Lambda}}{3M_{\text{Pl}}^2}, \qquad (C.5)$$

which is often expressed in terms of the Hubble function $H \equiv \dot{a}/a$. Due to the isotropic ansatz, the $ii$-components of the equations are all equivalent and therefore provide only one additional independent equation. The most common way to present it is by taking the traced Einstein equations and using (C.5) to arrive at the acceleration equation

$$\frac{\ddot{a}}{a} = -\frac{1}{6M_{\text{Pl}}^2}(\rho + \rho_{\Lambda} + 3p). \qquad (C.6)$$

In order to characterise different components of the cosmological fluid, one introduces the equation of state parameter $w$ defined as the ratio of pressure and energy density, $p = w\rho$. One can then rewrite the continuity equation (C.4) as

$$\frac{d\rho}{da} + 3(1+w)\frac{\rho}{a} = 0, \qquad (C.7)$$

which is solved by

$$\rho = \rho_0 a^{-3(1+w)}. \qquad (C.8)$$

Here, $\rho_0$ is a constant to be determined by boundary conditions. Inserting this solution into the Friedmann equation (C.5) with $k = 0$ and $\rho_{\Lambda} = 0$ gives the following evolution of the scale factor for $w \neq -1$
$$a \propto t^{\frac{2}{3(1+w)}}, \qquad (C.9)$$

The behaviour for $w = -1$ is

$$a \propto e^{Ht}, \qquad H \equiv \frac{\dot{a}}{a} = \text{const.} \qquad (C.10)$$

Moreover, from the acceleration equation (C.6) we infer that, in the absence of a bare cosmological constant, the expansion of the Universe accelerates for $w < -1/3$ and decelerates for $w > -1/3$. Three cases are of particular interest: $w = 0$ characterises non-relativistic matter, $w = 1/3$ corresponds to relativistic matter ('radiation'), and $w = -1$ describes a cosmological constant. The energy density for non-relativistic matter therefore scales as $a^{-3}$, i.e. inversely proportional to the volume. Non-relativistic matter experiences an additional redshift and its energy density goes as $a^{-4}$. None of these can explain an accelerated expansion of the Universe. In contrast, a cosmological constant which of course has constant energy density leads to an exponentially accelerated expansion.

It is common to measure the contributions of different fluid components in terms of their density parameters

$$\Omega_i = \frac{\rho_i}{3M_{\text{PL}}^2 H_0^2}, \qquad (C.11)$$

where $i$ stands for either radiation, non-relativistic matter or the cosmological constant and $H_0$ denotes the Hubble parameter at the present time. One also defines a curvature contribution, $\Omega_k = -\frac{k}{a^2 H_0^2}$, and in terms of these the Friedmann equation becomes

$$\Omega_{\text{rad}} + \Omega_{\text{mat}} + \Omega_{\Lambda} + \Omega_k = \frac{H^2}{H_0^2}. \qquad (C.12)$$

In particular, at the present time the density parameters add up to one. Latest observational data [31] suggest the following (approximate) values for the cosmological parameters at the present time

$$H_0 \sim 70 \text{ km s}^{-1} \text{ Mpc}^{-1}, \qquad \Omega_{\text{rad}} \sim 10^{-5}, \qquad \Omega_{\text{mat}} \sim 0.3, \\ \Omega_{\Lambda} \sim 0.7, \qquad \Omega_k < 10^{-3}. \qquad (C.13)$$

These values imply that our Universe is flat and dominated by a cosmological constant component. As already discussed above, the energy density corresponding to this cosmological constant is extremely small compared to energy scales of the Standard Model of Particle Physics. Since so far we lack an explanation for this small value, the curious energy component $\rho_{\Lambda}$ is often referred to as *dark energy*. On top of that, observations also show that the matter component $\Omega_{\text{m}}$ is not dominated by familiar baryonic matter, but rather mainly consists of an unknown *dark matter* component.

Consequently, under the assumption that GR is indeed the theory of gravity, we must accept that 95% of the Universe's energy content is not at all understood. Nevertheless, when this obscurity is ignored, the so-called $\Lambda$CDM model (GR with a cosmological constant and cold dark matter) fits the observational data very well.

## References

[1] Hinterbichler K 2012 *Rev. Mod. Phys.* **84** 671

[2] de Rham C 2014 *Living Rev. Relativ.* **17** 7

[3] Fierz M and Pauli W 1939 *Proc. R. Soc. A* **173** 211

[4] Pauli W and Fierz M 1939 *Helv. Phys. Acta* **12** 297
[5] Rosen N 1940 *Phys. Rev.* **57** 150

[6] Rosen N 1977 *Proc. Topics in Theoretical and Experimental Gravitation Physic* (Erice, New York, 1975) pp 273–94

[7] Aragone C and Deser S 1971 *Nuovo Cimento A* **3** 709

[8] Isham C J, Salam A and Strathdee J A 1971 *Phys. Rev.* D **3** 867

[9] van Dam H and Veltman M J G 1970 *Nucl. Phys.* B **22** 397

[10] Zakharov V I 1970 *J. Exp. Theor. Phys.* Lett. **12** 312

Zakharov V I 1970 *Pisma Zh. Eksp. Teor. Fiz.* **12** 447

[11] Iwasaki Y 1970 *Phys. Rev.* D **2** 2255

[12] Vainshtein A I 1972 *Phys. Lett. B* **39** 393

[13] Babichev E and Crisostomi M 2013 *Phys. Rev.* D **88** 084002

[14] Babichev E and Deffayet C 2013 *Class. Quantum Grav.* **30** 184001

[15] Boulware D G and Deser S 1972 *Phys. Rev.* D **6** 3368

[16] Zumino B 1970 *Lectures On Elementary Particles And Quantum Field Theory* vol 2 (Cambridge, MA: Brandeis University) pp 437–500

[17] Chamseddine A H, Salam A and Strathdee J A 1978 *Nucl. Phys.* B **136** 248

[18] Groot Nibbelink S, Peloso M and Sexton M 2007 *Eur. Phys. J. C* **51** 741

[19] Pitts J B 2016 *Ann. Phys.* **365** 73–90

[20] Maheshwari A 1972 *Nuovo Cimento A* **8** 319

[21] Ogievetsky V I and Polubarinov I V 1965 *Ann. Phys.* **35** 167–208

[22] Ogievetskii V I and Polubarinov I V 1966 *Sov. Phys. Dokl.* **11** 135–7

[23] Rubakov V A 2004 arXiv:hep-th/0407104

[24] Dubovsky S L 2004 *J. High Energy Phys.* JHEP10(2004)076

[25] Blas D, Comelli D, Nesti F and Pilo L 2009 *Phys. Rev.* D **80** 044025

[26] Comelli D, Nesti F and Pilo L 2013 *Phys. Rev.* D **87** 124021

[27] Comelli D, Nesti F and Pilo L 2013 *J. High Energy. Phys.* JHEP07(2013)161

[28] Comelli D, Nesti F and Pilo L 2014 *J. Cosmol. Astropart. Phys.* JCAP11(2014)018

[29] Riess A G et al (Supernova Search Team Collaboration) 1998 *Astron. J.* **116** 1009

[30] Perlmutter S et al (Supernova Cosmology Project Collaboration) 1999 *Astrophys. J.* **517** 565

[31] Adam R et al (Planck Collaboration) 2015 arXiv:1502.01582 [astro-ph.CO]

[32] Weinberg S 1989 *Rev. Mod. Phys.* **61** 1

[33] Zwicky F 1933 *Helv. Phys. Acta* **6** 110

[34] Rubin V C and Ford W K Jr 1970 *Astrophys. J.* **159** 379

[35] Arkani-Hamed N, Dimopoulos S and Dvali G R 1998 *Phys. Lett.* B **429** 263

[36] Antoniadis I, Arkani-Hamed N, Dimopoulos S and Dvali G R 1998 *Phys. Lett.* B **436** 257

[37] Randall L and Sundrum R 1999 *Nucl. Phys.* B **557** 79

[38] Dvali G R, Gabadadze G and Porrati M 2000 *Phys. Lett.* B **485** 208

[39] Dvali G R, Gabadadze G, Kolanovic M and Nitti F 2001 *Phys. Rev.* D **64** 084004

[40] Rubakov V A and Tinyakov P G 2008 *Phys. Usp.* **51** 759

[41] de Rham C, Hofmann S, Khoury J and Tolley A J 2008 *J. Cosmol. Astropart. Phys.* JCAP02 (2008)011

[42] de Rham C, Dvali G, Hofmann S, Khoury J, Pujolas O, Redi M and Tolley A J 2008 *Phys. Rev. Lett.* **100** 251603

[43] Dvali G R, Gabadadze G and Porrati M 2000 *Phys. Lett.* B **484** 112

[44] Gregory R, Rubakov V A and Sibiryakov S M 2000 *Phys. Rev. Lett.* **84** 5928

[45] Hassan S F, Hofmann S and von Strauss M 2011 *J. Cosmol. Astropart. Phys.* JCAP01(2011)020

[46] Hassan S F and Rosen R A 2011 *J High Energy Phys.* JHEP07(2011)009

[47] Arkani-Hamed N, Georgi H and Schwartz M D 2003 *Ann. Phys.* **305** 96

[48] Lee B W, Quigg C and Thacker H B 1977 *Phys. Rev.* D **16** 1519

[49] Chanowitz M S and Gaillard M K 1985 *Nucl. Phys.* B **261** 379

[50] Arkani-Hamed N, Cheng H-C, Luty M A and Mukohyama S 2004 *J. High Energy. Phys.* JHEP05(2004)074

[51] Creminelli P, Nicolis A, Papucci M and Trincherini E 2005 *J. High Energy. Phys.* JHEP09 (2005)003

[52] Deffayet C and Rombouts J W 2005 *Phys. Rev.* D **72** 044003

[53] Gabadadze G 2009 *Phys. Lett. B* **681** 89

[54] de Rham C 2010 *Phys. Lett. B* **688** 137

[55] de Rham C and Gabadadze G 2010 *Phys. Lett. B* **693** 334
[56] Hassan S F and Rosen R A 2011 *Phys. Lett. B* **702** 90

[57] de Rham C and Gabadadze G 2010 *Phys. Rev. D* **82** 044020

[58] de Rham C, Gabadadze G and Tolley A J 2011 *Phys. Rev. Lett.* **106** 231101

[59] Hassan S F and Rosen R A 2012 *Phys. Rev. Lett.* **108** 041101

[60] Arnowitt R L, Deser S and Misner C W 2008 *Gen. Relativ. Gravit.* **40** 1997

[61] Hassan S F, Rosen R A and Schmidt-May A 2012 *J. High Energy Phys. JHEP02(2012)026*

[62] Hassan S F and Rosen R A 2012 *J. High Energy Phys. JHEP04(2012)123*

[63] Hassan S F, Schmidt-May A and von Strauss M 2012 *Phys. Lett. B* **715** 335

[64] Golovnev A 2012 *Phys. Lett. B* **707** 404

[65] Kluson J 2012 *Phys. Rev. D* **86** 044024

[66] Kluson J 2012 arXiv:1209.3612 [hep-th]

[67] Alexandrov S, Krasnov K and Speziale S 2013 *J. High Energy Phys. JHEP06(2013)068*

[68] Alexandrov S 2014 *Gen. Relativ. Gravit.* **46** 1639

[69] Kugo T and Ohta N 2014 *Prog. Theor. Exp. Phys.* **2014** 043B04

[70] Damour T and Kogan I I 2002 *Phys. Rev. D* **66** 104024

[71] Damour T, Kogan I I and Papazoglou A 2002 *Phys. Rev. D* **66** 104025

[72] Blas D, Deffayet C and Garriga J 2006 *Class. Quantum Grav.* **23** 1697

[73] Berezhiani Z, Nesti F, Pilo L and Rossi N 2009 *J. High Energy Phys. JHEP07(2009)083*

[74] Berezhiani Z, Pilo L and Rossi N 2010 *Eur. Phys. J. C* **70** 305

[75] Deffayet C and Randjbar-Daemi S 2011 *Phys. Rev. D* **84** 044053

[76] Hassan S F and Rosen R A 2012 *J. High Energy Phys. JHEP02(2012)126*

[77] Butcher L M, Hobson M and Lasenby A 2009 *Phys. Rev. D* **80** 084014

[78] Deser S 2010 *Gen. Relativ. Gravit.* **42** 641

[79] Nieuwenhuizen P Van 1973 *Nucl. Phys. B* **60** 478

[80] Porrati M 2001 *Phys. Lett. B* **498** 92

[81] Higuchi A 1987 *Nucl. Phys. B* **282** 397

[82] Higuchi A 1989 *Nucl. Phys. B* **325** 745

[83] Kogan I I, Mouslopoulos S and Papazoglou A 2001 *Phys. Lett. B* **503** 173

[84] Deser S and Waldron A 2001 *Phys. Rev. Lett.* **87** 031601

[85] Deser S and Waldron A 2001 *Nucl. Phys. B* **607** 577

[86] Deser S and Waldron A 2001 *Phys. Lett. B* **508** 347

[87] Grisa L and Sorbo L 2010 *Phys. Lett. B* **686** 273

[88] Berkhahn F, Dietrich D D and Hofmann S 2010 *J. Cosmol. Astropart. Phys. JCAP11(2010)018*

[89] Berkhahn F, Dietrich D D and Hofmann S 2011 *J. Cosmol. Astropart. Phys. JCAP09(2011)024*

[90] Banados M, Gomberoff A, Rodrigues D C and Skordis C 2009 *Phys. Rev. D* **79** 063515

[91] Boulware D G and Deser S 1972 *Phys. Lett. B* **40** 227

[92] Date G 2010 arXiv:1010.2062 [gr-qc]

[93] DeWitt B S 1967 *Phys. Rev.* **160** 1113

[94] Hojman S A, Kuchar K and Teitelboim C 1976 *Ann. Phys.* **96** 88

[95] Volkov M S 2014 *Phys. Rev. D* **90** 124090

[96] Iglesias A and Kakushadze Z 2011 *Phys. Rev. D* **84** 084005

[97] Kakushadze Z 2013 *Universal J. Phys. Appl.* **7** 429

[98] Hassan S F and Kocic M private communication (to appear)

[99] Hassan S F, Kocic M and Schmidt-May A 2014 arXiv:1409.1909 [hep-th]

[100] Kluson J 2012 *J. High Energy Phys. JHEP01(2012)013*

[101] Deffayet C, Mourad J and Zahariade G 2013 *J. Cosmol. Astropart. Phys. JCAP01(2013)032*

[102] Deser S, Sandora M, Waldron A and Zahariade G 2014 *Phys. Rev. D* **90** 104043

[103] Buchbinder I L, Gitman D M, Krykhtin V A and Pershin V D 2000 *Nucl. Phys. B* **584** 615

[104] Bernard L, Deffayet C and von Strauss M 2015 *Phys. Rev. D* **91** 104013

[105] Bernard L, Deffayet C and von Strauss M 2015 *J. Cosmol. Astropart. Phys. JCAP038(2015)038*

[106] Guarato P and Durrer R 2014 *Phys. Rev. D* **89** 084016

[107] Cusin G, Durrer R, Guarato P and Motta M 2015 arXiv:1512.02131 [astro-ph.CO]

[108] Bernard L, Deffayet C, Schmidt-May A and von Strauss M 2015 arXiv:1512.03620 [hep-th]

[109] Lee B W, Quigg C and Thacker H B 1977 *Phys. Rev. Lett.* **38** 883

[110] Ostrogradski M 1850 *Mem. Ac. St. Petersburg VI* **4** 385

[111] Woodard R P 2015 *Scholarpedia* **10** 32242

[112] Nicolis A, Rattazzi R and Trincherini E 2009 *Phys. Rev. D* **79** 064036

[113] Mirbabayi M 2012 *Phys. Rev. D* **86** 084006
[114] de Rham C, Gabadadze G and Tolley A J 2012 *Phys. Lett. B* **711** 190

[115] Hassan S F, Schmidt-May A and von Strauss M 2013 *J. High Energy Phys. JHEP05(2013)086*

[116] Alberte L, Chamseddine A H and Mukhanov V 2011 *J. High Energy Phys. JHEP04(2011)004*

[117] Ondo N A and Tolley A J 2013 *J. High Energy Phys. JHEP11(2013)059*

[118] Chamseddine A H and Mukhanov V 2011 *J. High Energy Phys. JHEP08(2011)091*

[119] Kluson J 2012 *J. High Energy Phys. JHEP06(2012)170*

[120] Kluson J 2012 *Phys. Rev. D* **86** 124005

[121] de Rham C, Gabadadze G and Tolley A J 2011 *J. High Energy Phys. JHEP11(2011)093*

[122] Deser S and Waldron A 2013 *Phys. Rev. Lett.* **110** 111101

[123] Deser S, Izumi K, Ong Y C and Waldron A 2013 *Phys. Lett. B* **726** 544

[124] Deser S, Izumi K, Ong Y C and Waldron A 2013 arXiv:1312.1115 [hep-th]

[125] Deser S, Izumi K, Ong Y C and Waldron A 2015 *Mod. Phys. Lett. A* **30** 1540006

[126] Deser S, Waldron A and Zahariade G 2015 *Phys. Lett. B* **749** 144

[127] D'Amico G, de Rham C, Dubovsky S, Gabadadze G, Pirtskhalava D and Tolley A J 2011 *Phys. Rev. D* **84** 124046

[128] Gümrükçüoğlu A E, Lin C and Mukohyama S 2011 *J. Cosmol. Astropart. Phys* JCAP030 (2011)030

[129] Gümrükçüoğlu A E, Lin C and Mukohyama S 2012 *J. Cosmol. Astropart. Phys* JCAP03 (2012)006

[130] Vakili B and Khosravi N 2012 *Phys. Rev. D* **85** 083529

[131] De Felice A, Gümrükçüoğlu A E and Mukohyama S 2012 *Phys. Rev. Lett.* **109** 171101

[132] Fasiello M and Tolley A J 2012 *J. Cosmol. Astropart. Phys* JCAP11(2012)035

[133] De Felice A, Gümrükçüoğlu A E, Lin C and Mukohyama S 2013 *J. Cosmol. Astropart. Phys. JCAP05(2013)035*

[134] Fasiello M and Tolley A J 2013 *J. Cosmol. Astropart. Phys. JCAP12(2013)002*

[135] de Rham C, Matas A and Tolley A J 2014 *Class. Quantum Grav.* **31** 165004

[136] de Rham C, Matas A and Tolley A J 2015 *Class. Quantum Grav.* **32** 215027

[137] Matas A 2015 *Class. Quantum Grav.* **33** 075004

[138] Nojiri S and Odintsov S D 2012 *Phys. Lett. B* **716** 377

[139] Nojiri S, Odintsov S D and Shirai N 2013 *J. Cosmol. Astropart. Phys. JCAP05(2013)020*

[140] Kluson J, Nojiri S and Odintsov S D 2013 *Phys. Lett. B* **726** 918

[141] Bamba K, Makarenko A N, Myagky A N, Nojiri S and Odintsov S D 2014 *J. Cosmol. Astropart. Phys. JCAP01(2014)008*

[142] Nojiri S, Odintsov S D and Oikonomou V K 2015 *Phys. Rev. D* **92** 124059

[143] Hassan S F, Schmidt-May A and von Strauss M 2013 *Class. Quantum Grav.* **30** 184010

[144] Yamashita Y, De Felice A and Tanaka T 2014 *Int. J. Mod. Phys. D* **23** 3003

[145] de Rham C, Heisenberg L and Ribeiro R H 2015 *Class. Quantum Grav.* **32** 035022

[146] Baccetti V, Martin-Moruno P and Visser M 2013 *Class. Quantum Grav.* **30** 015004

[147] Hassan S F, Schmidt-May A and von Strauss M 2014 *Int. J. Mod. Phys. D* **23** 1443002

[148] Akrami Y, Hassan S F, Könnig F, Schmidt-May A and Solomon A R 2015 *Phys. Lett. B* **748** 37

[149] Aoki K and Maeda K 2014 *Phys. Rev. D* **90** 124089

[150] Bernard L and Blanchet L 2015 *Phys. Rev. D* **91** 103536

[151] Blanchet L and Heisenberg L 2015 *Phys. Rev. D* **91** 103518

[152] Blanchet L and Heisenberg L 2015 *J. Cosmol. Astropart. Phys. JCAP12(2015)026*

[153] Volkov M S 2012 *J. High Energy Phys. JHEP01(2012)035*

[154] von Strauss M, Schmidt-May A, Enander J, Mörtsell E and Hassan S F 2012 *J. Cosmol. Astropart. Phys. JCAP03(2012)042*

[155] Comelli D, Crisostomi M, Nesti F and Pilo L 2012 *J. High Energy Phys. JHEP03(2012)067*
Comelli D, Crisostomi M, Nesti F and Pilo L 2012 *J. High Energy Phys. JHEP06(2012)020*
(erratum)

[156] Park M and Sorbo L 2013 *Phys. Rev. D* **87** 024041

[157] Zhang Y L, Sasaki M and Yeom D h 2015 *J. High Energy Phys. JHEP04(2015)016*

[158] Akrami Y, Koivisto T S, Mota D F and Sandstad M 2013 *J. Cosmol. Astropart. Phys. JCAP10 (2013)046*

[159] Akrami Y, Koivisto T S and Solomon A R 2015 *Gen. Relativ. Gravit.* **47** 1838

[160] Melville S and Noller J 2016 *J. High Energy Phys. JHEP01(2016)094*

[161] Huang Q G, Ribeiro R H, Xing Y H, Zhang K C and Zhou S Y 2015 *Phys. Lett. B* **748** 356

[162] Heisenberg L 2015 *Phys. Rev. D* **92** 023525
[163] Cartan E 1933 C. R. Acad. Sci. Paris **196** 582-6

[164] Bekenstein J D 1993 Phys. Rev. D **48** 3641

[165] Schmidt-May A 2015 J. Cosmol. Astropart. Phys. JCAP01(2015)039

[166] de Rham C, Heisenberg L and Ribeiro R H 2014 Phys. Rev. D **90** 124042

[167] Solomon A R, Enander J, Akrami Y, Koivisto T S, Könnig F and Mörtsell E 2015 J. Cosmol. Astropart. Phys. JCAP04(2015)027

[168] Enander J, Solomon A R, Akrami Y and Mörtsell E 2015 J. Cosmol. Astropart. Phys. JCAP01 (2015)006

[169] Emir Gümçeküoğlu A, Heisenberg L and Mukohyama S 2015 J. Cosmol. Astropart. Phys. JCAP02(2015)022

[170] Heisenberg L 2015 Class. Quantum Grav. **32** 105011

[171] Gao X and Yoshida D 2015 Phys. Rev. D **92** 044057

[172] Gümçeküoğlu A E, Heisenberg L, Mukohyama S and Tanahashi N 2015 J. Cosmol. Astropart. Phys. JCAP04(2015)008

[173] Comelli D, Crisostomi M, Koyama K, Pilo L and Tasinato G 2015 J. Cosmol. Astropart. Phys. JCAP04(2015)026

[174] Lagos M and Noller J 2016 J. High Energy Phys. JHEP01(2016)023

[175] Deser S and Franklin J 2005 Am. J. Phys. **73** 261

[176] Volkov M S 2013 Class. Quantum Grav. **30** 184009

[177] Volkov M S 2015 Lecture Notes Phys. **892** 161

[178] Babichev E and Brito R 2015 Class. Quantum Grav. **32** 154001

[179] Gratia P, Hu W and Wyman M 2014 Phys. Rev. D **89** 027502

[180] Baccetti V, Martin-Moruno P and Visser M 2012 J. High Energy Phys. JHEP08(2012)148

[181] Sushkov S V and Volkov M S 2015 J. Cosmol. Astropart. Phys. JCAP06(2015)017

[182] Comelli D, Crisostomi M, Nesti F and Pilo L 2012 Phys. Rev. D **85** 024044

[183] Volkov M S 2012 Phys. Rev. D **85** 124043

[184] Enander J and Mörtsell E 2013 J. High Energy Phys. JHEP10(2013)031

[185] Brito R, Cardoso V and Pani P 2013 Phys. Rev. D **87** 124024

[186] Brito R, Cardoso V and Pani P 2013 Phys. Rev. D **88** 064006

[187] Brito R, Cardoso V and Pani P 2013 Phys. Rev. D **88** 023514

[188] Kodama H and Arraut I 2014 Prog. Theor. Exp. Phys. **2014** 023E02

[189] Babichev E and Fabbri A 2014 J. High Energy Phys. JHEP07(2014)016

[190] Babichev E and Fabbri A 2014 Phys. Rev. D **90** 084019

[191] Babichev E and Fabbri A 2014 Phys. Rev. D **89** 081502

[192] Enander J and Mortsell E 2015 J. Cosmol. Astropart. Phys. JCAP11(2015)023

[193] Ayón-Beato E, Higuita-Borja D and Méndez-Zavaleta J A 2016 Phys. Rev. D **93** 024049

[194] Deser S and Nepomechie R I 1984 Ann. Phys. **154** 396

[195] Babichev E and Fabbri A 2013 Class. Quantum Grav. **30** 152001

[196] Gregory R and Laflamme R 1993 Phys. Rev. Lett. **70** 2837

[197] Harmark T, Niarchos V and Obers N A 2007 Class. Quantum Grav. **24** R1

[198] 't Hooft G 1980 NATO Sci. Ser. B **59** 135

[199] Burgess C P 2013 arXiv:1309.4133 [hep-th]

[200] de Rham C, Gabadadze G, Heisenberg L and Pirtskhalava D 2011 Phys. Rev. D **83** 103516

[201] Akrami Y, Koivisto T S and Sandstad M 2013 J. High Energy Phys. JHEP03(2013)099

[202] Solomon A R 2015 arXiv:1508.06859 [gr-qc]

[203] Comelli D, Crisostomi M and Pilo L 2014 Phys. Rev. D **90** 084003

[204] De Felice A, Gümçeküoğlu A E, Mukohyama S, Tanahashi N and Tanaka T 2014 J. Cosmol. Astropart. Phys. JCAP06(2014)037

[205] Nersisyan H, Akrami Y and Amendola L 2015 Phys. Rev. D **92** 104034

[206] Akrami Y, Koivisto T S and Sandstad M 2013 arXiv:1302.5268 [astro-ph.CO]

[207] Könnig F, Patil A and Amendola L 2014 J. Cosmol. Astropart. Phys. JCAP03(2014)029

[208] Könnig F, Akrami Y, Amendola L, Motta M and Solomon A R 2014 Phys. Rev. D **90** 124014

[209] Cusin G, Durrer R, Guarato P and Motta M 2015 J. Cosmol. Astropart. Phys. JCAP05(2015)030

[210] Könnig F 2015 Phys. Rev. D **91** 104019

[211] Cusin G, Durrer R, Guarato P and Motta M 2015 J. Cosmol. Astropart. Phys. JCAP09(2015)043

[212] Comelli D, Crisostomi M and Pilo L 2012 J. High Energy Phys. JHEP06(2012)085

[213] Ratra B and Peebles P J E 1988 Phys. Rev. D **37** 3406

[214] de Rham C, Heisenberg L and Ribeiro R H 2013 Phys. Rev. D **88** 084058
[215] Berg M, Buchberger I, Enander J, Mörtsell E and Sjörs S 2012 *J. Cosmol. Astropart. Phys.* JCAP12(2012)021

[216] Solomon A R, Akrami Y and Koivisto T S 2014 *J. Cosmol. Astropart. Phys.* JCAP10(2014)066

[217] Lagos M and Ferreira P G 2014 *J. Cosmol. Astropart. Phys.* JCAP12(2014)026

[218] Enander J, Akrami Y, Mörtsell E, Renneby M and Solomon A R 2015 *Phys. Rev.* D **91** 084046

[219] Könnig F and Amendola L 2014 *Phys. Rev.* D **90** 044030

[220] Mörtsell E and Enander J 2015 *J. Cosmol. Astropart. Phys.* JCAP10(2015)044

[221] Deser S and Waldron A 2001 *Phys. Lett. B* **513** 137

[222] Deser S and Waldron A 2004 *Phys. Lett. B* **603** 30

[223] Deser S and Waldron A 2006 *Phys. Rev.* D **74** 084036

[224] Garidi T 2003 arXiv:hep-th/0309104

[225] Gabadadze G and Iglesias A 2008 *J. Cosmol. Astropart. Phys.* JCAP02(2008)014

[226] Deser S and Waldron A 2013 *Phys. Rev.* D **87** 087702

[227] Hinterbichler K 2015 *Phys. Rev.* D **91** 026008

[228] Zinoviev Y M 2007 *Nucl. Phys. B* **770** 83

[229] Joung E, Lopez L and Taronna M 2012 *J. High Energy Phys.* JHEP07(2012)041

[230] Joung E, Lopez L and Taronna M 2013 *J. High Energy Phys.* JHEP01(2013)168

[231] Deser S, Joung E and Waldron A 2013 *J. Phys. A: Math. Theor.* **46** 214019

[232] Deser S, Joung E and Waldron A 2012 *Phys. Rev.* D **86** 104004

[233] Deser S, Sandora M and Waldron A 2013 *Phys. Rev.* D **87** 101501

[234] Deser S, Sandora M and Waldron A 2013 *Phys. Rev.* D **88** 081501

[235] de Rham C, Hinterbichler K, Rosen R A and Tolley A J 2013 *Phys. Rev.* D **88** 024003

[236] Joung E, Li W and Taronna M 2014 *Phys. Rev. Lett.* **113** 091101

[237] Garcia-Saenz S and Rosen R A 2015 *J. High Energy Phys.* JHEP05(2015)042

[238] de Rham C and Renaux-Petel S 2013 *J. Cosmol. Astropart. Phys.* JCAP01(2013)035

[239] Hassan S F, Schmidt-May A and von Strauss M 2013 *Phys. Lett. B* **726** 834

[240] Hassan S F, Schmidt-May A and von Strauss M 2015 *Universe* **1** 92

[241] Stelle K S 1977 *Phys. Rev.* D **16** 953

[242] Bach R 1921 *Math. Z.* **9** 110

[243] Kaku M, Townsend P K and van Nieuwenhuizen P 1977 *Phys. Lett. B* **69** 304

[244] Fradkin E S and Tseytlin A A 1982 *Nucl. Phys. B* **201** 469

[245] Lee S C and van Nieuwenhuizen P 1982 *Phys. Rev.* D **26** 934

[246] Riegert R J 1984 *Phys. Lett. A* **105** 110

[247] Maldacena J 2011 arXiv:1105.5632 [hep-th]

[248] Mannheim P D 2012 *Found. Phys.* **42** 388

[249] Bender C M and Mannheim P D 2008 *Phys. Rev. Lett.* **100** 110402

[250] Bergshoeff E A, Hohm O and Townsend P K 2009 *Phys. Rev. Lett.* **102** 201301

[251] Stelle K S 1978 *Gen. Relativ. Gravit.* **9** 353

[252] Cusin G, Fumagalli J and Maggiore M 2014 *J. High Energy Phys.* JHEP09(2014)181

[253] Hassan S F, Schmidt-May A and von Strauss M 2015 *Class. Quantum Grav.* **33** 015011

[254] Garcia-Saenz S, Hinterbichler K, Joyce A, Mitsou E and Rosen R A 2016 *J. High Energy Phys.* JHEP02(2016)043

[255] D'Amico G, Gabadadze G, Hui L and Pirtskhalava D 2013 *Phys. Rev.* D **87** 064037

[256] De Felice A, Emir Gümrükçüoğlu A and Mukohyama S 2013 *Phys. Rev.* D **88** 124006

[257] Darabi F and Mousavi M 2015 *J. Cosmol. Astropart. Phys.* JCAP08(2015)009

[258] Boulanger N, Damour T, Gualtieri L and Henneaux M 2001 *Nucl. Phys. B* **597** 127

[259] Rham C d, Matas A, Ondo N and Tolley A J 2015 *Class. Quantum Grav.* **32** 175008

[260] Folkerts S, Pritzel A and Wintergerst N 2011 arXiv:1107.3157 [hep-th]

[261] Hinterbichler K 2013 *J. High Energy Phys.* JHEP10(2013)102

[262] Li W 2015 arXiv:1508.03246 [gr-qc]

[263] Lanczos C 1938 *Ann. Math.* **39** 842

[264] Lovelock D 1971 *J. Math. Phys.* **12** 498

[265] Lovelock D 1972 *J. Math. Phys.* **13** 874

[266] Deser S and Franklin J 2012 *Class. Quantum Grav.* **29** 072001

[267] Sisman T C, Gullu I and Tekin B 2012 *Phys. Rev.* D **86** 044041

[268] Paulos M F and Tolley A J 2012 *J. High Energy Phys.* JHEP09(2012)002

[269] Hinterbichler K and Rosen R A 2012 *J. High Energy Phys.* JHEP07(2012)047

[270] Hassan S F, Schmidt-May A and von Strauss M 2012 arXiv:1204.5202 [hep-th]
[271] Volkov M S 2012 *Phys. Rev. D* **86** 061502

[272] Zinoviev Y M 2013 *Class. Quantum Grav.* **30** 055005

[273] Kluson J 2014 *Eur. Phys. J. C* **74** 2985

[274] Bergshoeff E A, de Haan S, Hohm O, Merbis W and Townsend P K 2013 *Phys. Rev. Lett.* **111** 111102

Bergshoeff E A, de Haan S, Hohm O, Merbis W and Townsend P K 2013 *Phys. Rev. Lett.* **111** 259902

[275] Tamanini N, Saridakis E N and Koivisto T S 2014 *J. Cosmol. Astropart. Phys.* JCAP02 (2014)015

[276] de Rham C, Matas A and Tolley A J 2014 *Class. Quantum Grav.* **31** 025004

[277] Noller J, Scargill J H C and Ferreira P G 2014 *J. Cosmol. Astropart. Phys.* JCAP02(2014)007

[278] Noller J and Melville S 2015 *J. Cosmol. Astropart. Phys.* JCAP01(2015)003

[279] Soloviev V O 2014 arXiv:1410.0048 [hep-th]

[280] Scargill J H C, Noller J and Ferreira P G 2014 *J. High Energy Phys.* JHEP12(2014)160

[281] Noller J and Scargill J H C 2015 *J. High Energy Phys.* JHEP05(2015)034

[282] Hinterbichler K and Rosen R A 2015 *Phys. Rev. D* **92** 024030

[283] Nomura K and Soda J 2012 *Phys. Rev. D* **86** 084052

[284] Schmidt-May A and von Strauss M 2014 arXiv:1412.3812 [hep-th]

[285] Deser S and van Nieuwenhuizen P 1974 *Phys. Rev. D* **10** 411

[286] Deffayet C, Mourad J and Zahariade G 2013 *J. High Energy Phys.* JHEP03(2013)086

[287] Bañados M, Deffayet C and Pino M 2013 *Phys. Rev. D* **88** 124016

[288] de Rham C and Tolley A J 2015 *Phys. Rev. D* **92** 024024

[289] Porrati M 2002 *J. High Energy Phys.* JHEP04(2002)058

[290] Porrati M 2003 *Mod. Phys. Lett. A* **18** 1793

[291] Goon G, Hinterbichler K, Joyce A and Trodden M 2015 *J. High Energy Phys.* JHEP07(2015)101

[292] Bekaert X, Cnockaert S, Iazeolla C and Vasiliev M A 2005 arXiv:hep-th/0503128

[293] Taronna M 2012 arXiv:1209.5755 [hep-th]

[294] Didenko V E and Skvortsov E D 2014 arXiv:1401.2975 [hep-th]