# Massive mimetic cosmology

Adam R. Solomon, Valeri Vardanyan, Yashar Akrami

► To cite this version:

Adam R. Solomon, Valeri Vardanyan, Yashar Akrami. Massive mimetic cosmology. Physics Letters B, 2019, 794, pp.135-142. 10.1016/j.physletb.2019.05.045. hal-02066439

HAL Id: hal-02066439

https://hal.science/hal-02066439v1

Submitted on 26 Aug 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Massive mimetic cosmology

Adam R. Solomon<sup>a, *, </sup> Valeri Vardanyan<sup>b, c</sup> Yashar Akrami<sup>d, b</sup>

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

## ARTICLE INFO

Article history:

Received 14 March 2019

Received in revised form 24 May 2019

Accepted 28 May 2019

Available online 30 May 2019

Editor: H. Peiris

## ABSTRACT

We study the first cosmological implications of the mimetic theory of massive gravity recently proposed by Chamseddine and Mukhanov. This is a novel theory of ghost-free massive gravity which additionally contains a mimetic dark matter component. In an echo of other modified gravity theories, there are self-accelerating solutions which contain a ghost instability. In the ghost-free region of parameter space, the effect of the graviton mass on the cosmic expansion history amounts to an effective negative cosmological constant, a radiation component, and a negative curvature term. This allows us to place constraints on the model parameters—the graviton mass and the Stückelberg vacuum expectation value—by insisting that the effective radiation and curvature terms be within observational bounds. The late-time acceleration must be accounted for by a separate positive cosmological constant or other dark energy sector. We impose further constraints at the level of perturbations by demanding linear stability. We comment on the possibility of distinguishing this theory from ACDM with current and future large-scale structure surveys.

© 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ([http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)). Funded by SCOAP³.

## 1. Introduction

Chamseddine and Mukhanov have recently proposed [1,2] a novel ghost-free theory of massive gravity in which one of the four Stückelberg scalars is constrained in the same way as in the mimetic theory of dark matter [3], spontaneously breaking Lorentz invariance. In this Letter, we study the immediate implications of this mimetic massive gravity for cosmological theory and observation.

From a field-theoretic perspective, general relativity is the unique theory (in four spacetime dimensions) of a massless spin-2 particle, or graviton. It is therefore natural to ask whether it is possible to endow the graviton with a non-zero mass, and what sort of theoretical structures would result [4]. A closely related line of inquiry asks whether it is possible for two or more gravitons to interact [5]. Most nonlinear realizations of such theories suffer from the so-called Boulware-Deser ghost instability [6]. The past decade has seen the construction of models which avoid this instability, allowing for the construction of ghost-free theories of

massive gravity [7–13] and bimetric and multimetric gravity [12, 14, 15]. We refer the reader to the reviews [16,17] on massive gravity and [18,19] on bimetric gravity. The theory of mimetic massive gravity proposed in Refs. [1,2] takes a new and alternative path to a ghost-free nonlinear theory of massive gravity.

A generic theory of massive gravity propagates six degrees of freedom, which should be thought of as the five helicity states of a massive graviton plus an additional, ghostly scalar. The easiest way to understand the degrees-of-freedom counting is to observe that a graviton mass breaks diffeomorphism invariance. This is a gauge symmetry and so can be restored by the addition of four Stückelberg scalars $\Phi^A$, which propagate in addition to the two (now potentially massive) tensor modes of general relativity.

As an illustration, consider a Lorentz-invariant theory of massive gravity. In order to construct non-trivial, non-derivative interactions for the metric, one requires a second “reference” metric. The simplest choice for this metric is that of flat space, $\eta_{\mu\nu}$, but the addition of this prior geometry breaks diffeomorphism invariance; for instance, there are preferred coordinate systems in which $\eta_{\mu\nu} = \text{diag}(-1, 1, 1, 1)$. But diffeomorphism invariance is simply a redundancy in description, and can be restored by the addition of redundant variables, i.e., replacing $\eta_{\mu\nu} \rightarrow \eta_{AB}\partial_{\mu}\Phi^{A}\partial_{\nu}\Phi^{B}$, where $\eta_{AB} = \text{diag}(-1, 1, 1, 1)$ and the four fields $\Phi^A$ transform as

<footnote><sup>a</sup> Department of Physics & McWilliams Center for Cosmology, Carnegie Mellon University, Pittsburgh, PA 15213, USA</footnote>

<footnote><sup>b</sup> Lorentz Institute for Theoretical Physics, Leiden University, P.O. Box 9506, 2300 RA Leiden, the Netherlands</footnote>

<footnote><sup>c</sup> Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, the Netherlands</footnote>

<footnote><sup>d</sup> Département de Physique, École Normale Supérieure, PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: adamsolo@andrew.cmu.edu (A.R. Solomon),
vardanyan@lorentz.leidenuniv.nl (V. Vardanyan), akrami@ens.fr (Y. Akrami).</footnote>
spacetime scalars. One can always, by means of a diffeomorphism, choose the unitary gauge in which $\Phi^A = x^A$, and we recover the original description of the theory in terms of a symmetry-breaking reference metric. Generic interaction terms for the graviton, e.g., generic functions of $g^{\mu\alpha}\eta_{AB}\partial_\alpha\Phi^A\partial_\nu\Phi^B$, will lead to dynamics for each of these four scalars, in addition to the two modes of general relativity, for a total of six degrees of freedom.

At the linear level, i.e., linearizing the metric about flat space in unitary gauge, $g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}$ and $\Phi^A = x^A$, we find that one of the six degrees of freedom leads to a ghost instability unless we specifically arrange the mass term into the Fierz-Pauli form, $\mathcal{L}_{\text{mass}} \sim h_{\mu\nu}^2 - h^2$, in which case the dynamics of the ghostly mode take the form of a total derivative. Continuing this procedure at higher orders in perturbation theory—i.e., continually packaging ghostly operators into total derivative structures—leads uniquely to the non-linear massive gravity theory of de Rham, Gabadadze, and Tolley (dRGT) [8,9].

The recent proposal of Chamseddine and Mukhanov takes a novel alternative approach to eliminating the dangerous ghostly mode [1,2]. Noticing that the ghost can be associated to the $\Phi^0$ Stückelberg mode, they propose imposing the constraint $g^{\mu\nu}\partial_\mu\Phi^0\partial_\nu\Phi^0 = -1$. This is motivated by a similar construction known as mimetic gravity [3], in which the constrained scalar winds up behaving like dark matter.<fnref n="1" /> Mimetic massive gravity takes this constrained scalar to be one of the Stückelberg modes of a massive graviton, eliminating the ghost. They propose the following action, designed to ensure stability at the linear level (notice that the mass term is not of the Fierz-Pauli form),

$$S = \int d^4x \sqrt{-g} \left[ \frac{M_{\text{Pl}}^2}{2} R + \frac{m^2 M_{\text{Pl}}^2}{8} \left( \frac{1}{2}\bar{h}^2 - \bar{h}_{AB}^2 \right) + \lambda(X+1) \right] + S_{\text{matter}} \quad (1)$$

with $X \equiv g^{\mu\nu}\partial_\mu\Phi^0\partial_\nu\Phi^0$, and

$$\bar{h}^{AB} \equiv g^{\mu\nu} \partial_{\mu} \Phi^{A} \partial_{\nu} \Phi^{B} - \eta^{AB}. \quad (2)$$

Internal indices (given by capital Roman letters) are raised and lowered with the Minkowski metric. The field equations are<fnref n="2" />

$$G_{\mu\nu} = \frac{1}{M_{\text{Pl}}^2} T_{\mu\nu} - \frac{2\lambda}{M_{\text{Pl}}^2} \partial_{\mu} \Phi^0 \partial_{\nu} \Phi^0 + \frac{m^2}{2} \left( \bar{h}_{AB} - \frac{1}{2}\bar{h}\eta_{AB} \right) \left( \partial_{\mu} \Phi^A \partial_{\nu} \Phi^B - \frac{1}{4}\bar{h}^{AB} g_{\mu\nu} \right), \quad (3)$$

$$0 = \nabla_{\mu} \left[ \frac{2\lambda}{M_{\text{Pl}}^2} \partial^{\mu} \Phi^0 \delta_A^0 - \frac{m^2}{2} \left( \bar{h}_{AB} - \frac{1}{2}\bar{h}\eta_{AB} \right) \partial^{\mu} \Phi^B \right], \quad (4)$$

$$X = -1. \quad (5)$$

The last of these aligns $\Phi^0$ with the lapse of $g_{\mu\nu}$. An upshot of this construction is that the constrained mode behaves as a pressureless fluid, i.e., this theory provides a (mimetic) dark matter candidate [1,2].<fnref n="3" />

We end this section by making a connection with the existing literature on Lorentz-violating massive gravity and demonstrating the absence of certain well-known features of Lorentz-invariant massive gravity, namely the van Dam-Veltman-Zakharov (vDVZ) discontinuity [22,23] and the Higuchi bound [24]. The vDVZ discontinuity refers to the failure of linearized Lorentz-invariant massive gravity to reduce to general relativity in the massless limit; this requires nonlinear effects in order to restore general relativity in the Newtonian limit [25,26]. The Higuchi bound is a stability bound for massive gravity on de Sitter space, placing a lower bound on the graviton mass, $m^2 \ge 2H^2$, with H the Hubble rate. It is well-known that breaking Lorentz invariance changes both of these conclusions dramatically [27,28].

At the level of linear perturbations around flat space, the general SO(3)-invariant mass term in unitary gauge ($\Phi^A = x^A$) can be written as [27]

$$\mathcal{L}_{\text{mass}} = \frac{1}{8} M_{\text{Pl}}^2 \left( m_0^2 h_{00}^2 + 2m_1^2 h_{0i}^2 - m_2^2 h_{ij}^2 + m_3^2 h_{ii}^2 - 2m_4^2 h_{00} h_{ii} \right). \quad (6)$$

The linearized mass term in eq. (1) in unitary gauge is (treating $\lambda$ as first-order)

$$\mathcal{L}_{\text{mass}} = \frac{m^2 M_{\text{Pl}}^2}{8} \left( -\frac{1}{2}h_{00}^2 + 2h_{0i}^2 - h_{ij}^2 + \frac{1}{2}h_{ii}^2 - h_{00}h_{ii} \right) + \lambda h_{00}. \quad (7)$$

The $\lambda$ equation of motion sets $h_{00} = 0$, which we can impose in the action⁴ to find

$$m_0^2 = m_4^2 = 0, \quad m_1^2 = m_2^2 = 2m_3^2 = 1. \quad (8)$$

This allows us to easily make contact with the existing literature on Lorentz-violating massive gravity. The analysis of Ref. [27] shows that for these $m_i$ parameters, the Newtonian limit is the usual one, while the vDVZ discontinuity is absent. The analogue of the Higuchi bound in Lorentz-violating massive gravity was derived in Ref. [28], and for our values of the $m_i$ parameters, it reduces simply to $H^2 > 0$, which is trivially satisfied.

## 2. Flat-space perturbations

In this section, we briefly review the behavior of perturbations about flat space in mimetic massive gravity, as discussed in Refs. [1,2]. This will place stability conditions on the theory which will be relevant when we move to cosmological solutions.

The equations of motion (3)–(5) in vacuum are solved by⁵

$$g_{\mu\nu} = \eta_{\mu\nu}, \quad \Phi^A = x^A, \quad \lambda = 0. \quad (9)$$

ple, since the mimetic dark matter only interacts gravitationally with the Standard Model, we do not expect to have a thermal production mechanism, in contrast to many traditional dark matter scenarios such as WIMPs. Indeed, when the theory is shift-symmetric in $\Phi^0$, the energy density of this component is set entirely by an integration constant and so is determined by initial conditions. It may also be necessary to tune the parameters of the model in order to obtain the right values of the dark matter density over the entire cosmic history, and higher-derivative effective field theory corrections play an important role [21]. We refer the reader to, e.g., Ref. [21] for discussions of the constraints that early-Universe considerations place on the properties and evolution of mimetic dark matter throughout cosmic history.

⁴ This is justified because, on shell, the $h_{00}$ equation of motion simply sets the value of $\lambda$, while $h_{00}$ drops out of the $h_{ij}$ equations of motion. The dynamics are therefore equivalent.

⁵ This is the only solution that is manifestly invariant under rotations, i.e., with $g_{\mu\nu} = \text{diag}(-1, 1, 1, 1)$ and $\Phi^A = \{\varphi(t), \beta x^i\}$. A priori it may be possible to have flat solutions with inhomogeneous Stückelbergs $\Phi^A$, or equivalently solutions with $\Phi^A = x^A$ and $g_{\mu\nu} = \eta_{\mu\nu}$ with $\eta_{\mu\nu}$ written in a nonstandard coordinate system, but we do not consider these here.

<footnote n="1">¹ For an earlier construction in which a constrained scalar mimics dark matter and dark energy, and which contains mimetic dark matter as a subset, see Ref. [20].</footnote>

<footnote n="2">² Note the sign differences between the right-hand side of the Einstein equations and the corresponding equation in Ref. [1], which is due to the mostly positive metric convention we employ.</footnote>

<footnote n="3">³ One should note that the phenomenology of mimetic dark matter is still in the early stages of development compared to traditional particle dark matter models such as weakly interacting massive particles (WIMPs) or axions, and it is premature to consider mimetic gravity as a serious alternative to those models. For example,</footnote>
We expand the action (1) to quadratic order around the Minkowski solution, focusing on scalar modes,

$$g_{00} = -(1 + 2\phi), \quad (10)$$

$$g_{0i} = \partial_i B, \quad (11)$$

$$g_{ij} = (1 - 2\psi)\delta_{ij} + 2\partial_i \partial_j E, \quad (12)$$

$$\Phi^A = x^A + \{\pi^0, \partial^i \pi\}, \quad (13)$$

$$\lambda = \delta\lambda. \quad (14)$$

Three of these fields—$\phi$, $B$, and $\delta\lambda$—are auxiliary, as they appear without time derivatives in the action, and so can be integrated out using their equations of motion. Note that the auxiliary structure is precisely the same as in general relativity, since the mass term and Lagrange multiplier do not introduce any derivatives of the metric.

We can use diffeomorphism invariance to remove a further two modes. When gauge fixing at the level of the action, one must take care to only eliminate variables whose equations of motion are contained in the equations of motion of the remaining variables, otherwise we will lose information after picking a gauge. Following the procedure of Refs. [29,30], we see that we can safely take $\pi^0$ and one of $(E, \pi)$ to vanish. Picking unitary gauge, $\pi^0 = \pi = 0$, we obtain the flat-space quadratic action (in Fourier space),

$$\delta_2 S = \int dt M_{\mathrm{Pl}}^2 \left( -\vec{\dot{X}}^\mathrm{T} \mathbb{K} \vec{\dot{X}} + \vec{X}^\mathrm{T} \left( k^2 \mathbb{G} + m^2 \mathbb{M} \right) \vec{X} \right), \quad (15)$$

where $\vec{\dot{X}} \equiv (\psi, k^2 E)$ and the matrices $\mathbb{K}$, $\mathbb{G}$, and $\mathbb{M}$ are given by

$$\mathbb{K} = \begin{pmatrix} 3 + \frac{4k^2}{m^2} & 1 \\ 1 & 0 \end{pmatrix}, \quad (16)$$

$$\mathbb{G} = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}, \quad (17)$$

$$\mathbb{M} = \frac{1}{4} \begin{pmatrix} 3 & 1 \\ 1 & -1 \end{pmatrix}. \quad (18)$$

As described in Ref. [2], this system can be diagonalized by replacing $\psi$ with the Lagrange multiplier $\delta\lambda$, which we had previously integrated out using

$$\delta\lambda = \frac{M_{\mathrm{Pl}}^2}{4} \left[ (4k^2 + 3m^2)\psi + k^2 m^2 E \right], \quad (19)$$

to find

$$\delta_2 S = \int dt \frac{1}{4k^2 + 3m^2} \left[ k^4 m^2 M_{\mathrm{Pl}}^2 (\dot{E}^2 - (k^2 + m^2)E^2) - \frac{1}{M_{\mathrm{Pl}}^2} \left( \frac{16}{m^2} \delta\dot{\lambda}^2 - 4\delta\lambda^2 \right) \right]. \quad (20)$$

If we take $m^2 > 0$, we can canonically normalize,

$$\delta\lambda_c \equiv \frac{4}{m M_{\mathrm{Pl}} \sqrt{2k^2 + \frac{3}{2}m^2}} \delta\lambda, \quad (21)$$

$$E_c \equiv \frac{m M_{\mathrm{Pl}} k^2}{\sqrt{2k^2 + \frac{3}{2}m^2}} E, \quad (22)$$

to obtain the final action,

$$\delta_2 S = \int dt \left[ \frac{1}{2} \dot{E}_c^2 - \frac{1}{2} (k^2 + m^2) E_c^2 - \frac{1}{2} \delta\lambda_c^2 + \frac{1}{8} m^2 \delta\lambda_c^2 \right]. \quad (23)$$

The only dynamical degree of freedom here is $E_c$, which is healthy and has mass $m$. The field $\delta\lambda_c$ has the wrong sign on both

its kinetic and mass terms, but does not propagate due to the absence of a gradient term; its equation of motion,

$$\ddot{\delta\lambda}_c + \frac{m^2}{4} \delta\lambda_c = 0, \quad (24)$$

leads to a dispersion relation $\omega^2 = m^2/4$ and is solved simply by [2]

$$\delta\lambda_c = C(\vec{x}) \sin\left(\frac{mt}{2}\right) + D(\vec{x}) \cos\left(\frac{mt}{2}\right), \quad (25)$$

where $C$ and $D$ are space-dependent constants of integration. The authors of Ref. [2] identify this mode with the mimetic dark matter.<fnref n="6" />

When we discuss cosmology in the next section, we will find ourselves tempted by the possibility of taking $m^2 < 0$. A priori this is merely a parameter choice, but the flat-space analysis shows why this would be a poor decision. By looking at the action (20), we see that, for negative $m^2$, the overall sign in front of the action flips depending on whether $k^2 > 3|m^2|/4$ or $k^2 < 3|m^2|/4$, a sign of pathological behavior. In particular, for scales $k^2 > 3|m^2|/4$, upon canonically normalizing we find the action (23) with an overall minus sign, so that the dynamical mode $E_c$ is a ghost.

### 3. Cosmological solutions

In this section we investigate Friedmann-Lemaître-Robertson-Walker (FLRW) cosmological solutions of mimetic massive gravity. Consider the homogeneous and isotropic ansatz

$$g_{\mu\nu} = \mathrm{diag}(-1, a(t)^2 \delta_{ij}), \quad (26)$$

$$\Phi^A = \{\varphi(t), \beta x^i\}. \quad (27)$$

In principle one could allow $\beta$ to depend on time, but this breaks homogeneity and isotropy as it induces $\vec{x}$-dependent terms in the stress-energy tensor of the Stückelberg fields. Note that on-shell, the Lagrange multiplier enforces $\varphi = t$ (up to a constant). We will include a general matter sector with density $\rho$ and pressure $p$. We will find this sector needs to contain a cosmological constant, much like in general relativity, but does not need to include dark matter, as this role can be played by the mimetic dark matter (which is an exactly pressureless perfect fluid).

The Einstein and scalar equations of motion are

$$3H^2 = \frac{\rho}{M_{\mathrm{Pl}}^2} - \frac{2\lambda}{M_{\mathrm{Pl}}^2} - \frac{3m^2}{16} \left( \frac{\beta^4}{a^4} - 6\frac{\beta^2}{a^2} + 5 \right), \quad (28)$$

$$2\dot{H} + 3H^2 = -\frac{p}{M_{\mathrm{Pl}}^2} - \frac{m^2}{16} \left( 3 - \frac{\beta^4}{a^4} - 2\frac{\beta^2}{a^2} \right), \quad (29)$$

$$0 = \frac{d}{dt} \left\{ a^3 \left[ \frac{3m^2}{4} \left( 1 - \frac{\beta^2}{a^2} \right) + \frac{2\lambda}{M_{\mathrm{Pl}}^2} \right] \right\}. \quad (30)$$

We can solve for $\lambda$ by integrating the $\Phi^0$ equation of motion (30), finding

<footnote n="6">⁶ See Ref. [2] for an argument for why this mode is not a ghost, despite having an overall wrong-sign action. In principle, one might worry that when quantizing or considering nonlinearities, a coupling will be induced between $\delta\lambda_c$ and other fields which will lead to an Ostrogradski instability. On the other hand, due to the lack of a gradient term this mode is not a propagating degree of freedom in the usual sense. We will remain agnostic about this question and limit ourselves to considerations of classical, linear stability, which this system clearly satisfies for $m^2 > 0$. See, e.g., Refs. [31,32] for detailed discussions of classical and quantum properties of modes lacking a gradient term.</footnote>
$$ -\frac{2\lambda}{M_{\text{Pl}}^2} = \frac{C}{a^3} + \frac{3m^2}{4} \left(1 - \frac{\beta^2}{a^2}\right), \qquad (31) $$

where C is an integration constant. Plugging this into the Friedmann equation (28), we obtain

$$ 3H^2 = \frac{\rho}{M_{\text{Pl}}^2} + \frac{C}{a^3} - \frac{3m^2}{16} \left(1 - \frac{\beta^2}{a^2}\right)^2. \qquad (32) $$

Note that the contribution from $\lambda$ exactly cancels out that from the last term of the Einstein equation (3), so the very simple form for $\rho_{\varphi} \equiv -3m^2 M_{\text{Pl}}^2 (1 - \beta^2/a^2)^2/16$ is entirely due to the term proportional to $g_{\mu\nu}$ in the stress tensor. The integration constant provides a dust-like contribution to the Friedmann equation, which is to be expected as this is a theory of mimetic dark matter.

We can get a better sense of the physical picture by expanding out the Friedmann equation and absorbing the mimetic dark matter C into $\rho$, finding

$$ 3H^2 = \frac{\rho}{M_{\text{Pl}}^2} - \frac{3m^2}{16} \left( \frac{\beta^4}{a^4} - \frac{2\beta^2}{a^2} + 1 \right). \qquad (33) $$

For $m^2 > 0$ ($m^2 < 0$), we see that the mass term generates an effective negative (positive) cosmological constant, an effective negative (positive) curvature, and an effective radiation component with negative (positive) energy density. Note that these add on to any cosmological constant, radiation, and curvature already present cosmologically; for example, while we have assumed a flat cosmology as our ansatz, observational bounds on spatial curvature will constrain the sum of any pre-existing curvature and the curvature-like term generated by the graviton mass.

Note that for $m^2 < 0$ we have late-time acceleration, with $\Lambda_{\text{eff}} = 3|m^2|/16$. However, as discussed in the previous section, we need $m^2 > 0$ in order to avoid a ghost around flat space. This is reminiscent of the situation in the Dvali-Gabadadze-Porrati (DGP) model [33], where one branch of solutions has self-accelerating cosmological expansion [34,35] but is plagued by a ghost [36,37], while the other branch is healthy but cannot account for cosmic acceleration.

Let us assume that the energy density $\rho$ in eq. (33) contains dust (including the mimetic dark matter), radiation, and dark energy components. Then, in terms of the density parameters,

$$ \Omega_{i,0} = \frac{\rho_{i,0}}{3M_{\text{Pl}}^2 H_0^2}, \qquad (34) $$

the components of the Friedmann equation which are modified by mimetic massive gravity are

$$ \Omega_{\Lambda,0} = \bar{\Omega}_{\Lambda,0} - \frac{m^2}{16H_0^2} \qquad (35) $$

$$ \Omega_{K,0} = \frac{m^2}{8H_0^2} \beta^2, \qquad (36) $$

$$ \Omega_{r,0} = \bar{\Omega}_{r,0} - \frac{m^2}{16H_0^2} \beta^4, \qquad (37) $$

where $\bar{\Omega}_{\Lambda,0}$ and $\bar{\Omega}_{r,0}$ are the densities associated to dark energy and Standard Model radiation. Using observational bounds on the curvature and radiation densities, we can place constraints on the model parameters $m^2$ and $\beta$. We will not consider any bounds coming from the presence of the effective cosmological constant, even though it contributes a negative and potentially large (if $m^2 \gg H_0$) amount to $\Omega_{\Lambda,0}$. Particle physics also predicts a large (and potentially negative) vacuum energy, and since we are not worrying about that, it seems inconsistent to worry about the

contribution from mimetic massive gravity. One might expect that whatever solves the former problem will also solve the latter.<fnref n="7" />

We will use observational constraints on $\Omega_{K,0}$ and $\Omega_{r,0}$ to bound our two free parameters, $m^2$ and $\beta$. Planck 2018 constrains $\Omega_{K,0} = 0.0007 \pm 0.0019$, which we parametrize as $|\Omega_{K,0}| < \delta_K$, with $\delta_K \sim 0.003$ [39]. We will take this to be a constraint on the contribution from mimetic massive gravity alone,

$$ \frac{m^2}{8H_0^2} \beta^2 < \delta_K. \qquad (38) $$

We remind the reader that what we are really bounding is the sum of the mimetic massive gravity contribution and any “bare” curvature, but unless there is significant tuning between these two, we can simply take this as a constraint on the mimetic massive gravity piece alone.

To bound the mimetic contribution to the radiation density, we will use constraints from big bang nucleosynthesis (BBN). At the time of BBN, radiation dominates. The exact value of the Hubble rate at the time of nucleosynthesis, which depends on the radiation density, determines the freeze-out abundance of neutrons and therefore the primordial abundance of helium-4, which is subject to tight observational bounds. The constraints are conveniently phrased in terms of the “speed-up factor” $\zeta \equiv H/\dot{H}$, where $H$ and $\dot{H}$ are the Hubble rate and its expected value, respectively, at the time of BBN. The difference between the observed and predicted helium-4 abundance, $|\Delta Y_P|$, is related to the speed-up factor by [40]

$$ \Delta Y_P = 0.08(\zeta^2 - 1). \qquad (39) $$

Current observational bounds imply [41]

$$ |\Delta Y_P| \lesssim 0.01. \qquad (40) $$

Comparing the Friedmann equation (33) with and without the mimetic radiation contribution, and focusing on radiation domination, we find

$$ \zeta^2 - 1 = - \frac{m^2 \beta^4}{16 \bar{\Omega}_{r,0} H_0^2}, \qquad (41) $$

where the value for the present-day radiation density associated to photons and neutrinos, $\bar{\Omega}_{r,0} \sim 10^{-4}$, is determined entirely by the CMB temperature and the effective number of neutrino species and is therefore not dependent on our modification of gravity.<fnref n="8" /> Combining this with eq. (40) we arrive at the constraint

$$ \frac{m^2}{16H_0^2} \beta^4 < \delta_r, \qquad (42) $$

where

$$ \delta_r \equiv \frac{\max(|\Delta Y_P|)\bar{\Omega}_{r,0}}{0.08} \approx \mathcal{O}(10^{-5}). \qquad (43) $$

We can rewrite our constraints (38) and (42) as inequalities for $m/H_0$ and $\beta$ alone in two different régimes,

$$ \frac{m}{H_0} < \begin{cases} \frac{\sqrt{8\delta_K}}{\beta}, & \beta < \sqrt{\frac{2\delta_r}{\delta_K}} \\ \frac{4\sqrt{\delta_r}}{\beta^2}, & \beta > \sqrt{\frac{2\delta_r}{\delta_K}} \end{cases} \qquad (44) $$

These are plotted in Fig. 1.

<footnote n="7">⁷ See Ref. [38] for a proposed solution to the cosmological constant problem in the context of Lorentz-violating massive gravity, which is closely related to mimetic massive gravity.</footnote>

<footnote n="8">⁸ See Ref. [42] for a measurement of the CMB temperature.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Upper limits on $m/H_0$ and $\beta$ for $(\delta_K, \delta_r) = (0.003, 10^{-5})$.</figcaption></figure>

Finally, we note that the strong-coupling scale for this theory is of order $\Lambda_2 = \sqrt{m}M_{Pl}$ [2]. If $m$ is of order the present-day Hubble scale, $m \sim 10^{-33}$ eV, then the strong coupling scale is $\Lambda_2 \sim$ meV, i.e., the theory breaks down slightly below the millimeter scale. As we see from eq. (44), for sufficiently small $\beta$, $m$ could potentially be much larger than $H_0$, leading to a correspondingly larger strong-coupling scale.

## 4. Cosmological perturbations

As we have seen, at the background level, cosmological solutions in mimetic massive gravity do not differ appreciably from $\Lambda$CDM. We therefore proceed to study cosmological perturbations around the FLRW background. This will tell us how cosmological large-scale structure (LSS) evolves in this theory in comparison to $\Lambda$CDM. Since mimetic massive gravity differs from general relativity, we would expect modifications to the gravitational Poisson equation and the slip relation, which could in principle allow for observational tests of this alternative model against $\Lambda$CDM and distinguish the two using the current and future LSS surveys. However, as we will see, stability of cosmological perturbations and the bounds (44) place strong constraints on the model which suggest that this theory should be observationally indistinguishable from GR in the linear régime.

### 4.1. Stability bound

We begin by studying the stability of cosmological perturbations using the second-order action formalism. Since, as discussed in section 3, this theory does not possess ghost-free self-accelerating solutions, we include a cosmological constant, although it will not affect any of the results in this section. Since the theory already contains a pressureless fluid, namely the mimetic dark matter, we need not introduce an additional matter field. Our analysis is therefore valid for all times after matter-radiation equality.

We define the linearized metric, Stückelberg fields, and La-grange multiplier as

$$ds^2 = -(1 + 2\phi)dt^2 + 2a\partial_i B dt dx^i + a^2 [(1 - 2\psi)\delta_{ij} + 2\partial_i \partial_j E] dx^i dx^j, \quad (45)$$

$$\Phi^0 = t + \pi^0, \quad (46)$$

$$\Phi^i = \beta (\chi^i + \partial^i \pi), \quad (47)$$

$$\lambda = \bar{\lambda} + \delta\lambda, \quad (48)$$

where we are restricting ourselves to scalar perturbations, and $\bar{\lambda}$ is the background value given in eq. (31). The calculation of the

quadratic action proceeds analogously to the flat-space case discussed in section 2. Expanding the action (1) (with a cosmological constant) to quadratic order in perturbations, we find that the variables $\phi$, $B$, and $\delta\lambda$ are auxiliary—that is, they appear without time derivatives—and can therefore be integrated out using their equations of motion. To safely fix a gauge at the level of the action, we again follow the procedure of Refs. [29,30], finding that we can eliminate one each of $(\psi, \pi^0)$ and $(E, \pi)$. We will choose to work in unitary gauge, $\pi^0 = \pi = 0$, so that $\Phi^A = (t, \beta x^i)$ is unperturbed. The final action, in Fourier space and after integrations by parts, is

$$\delta_2 S = \int dt M_{Pl}^2 a^3 \left( -\dot{\vec{\chi}}^T \mathbb{K} \dot{\vec{\chi}} + \vec{\chi}^T \left( \frac{k^2}{a^2} \mathbb{G} + m^2 \mathbb{M} \right) \vec{\chi} \right), \quad (49)$$

where $\vec{\chi} \equiv (\psi, k^2 E)$ and the matrices $\mathbb{K}$, $\mathbb{G}$, and $\mathbb{M}$ are given by

$$\mathbb{K} = \begin{pmatrix} 3 - \frac{8a^2}{\beta^2 - 3a^2} \frac{k^2}{m^2 \beta^2} & 1 \\ 1 & 0 \end{pmatrix} \quad (50)$$

$$\mathbb{G} = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \quad (51)$$

$$\mathbb{M} = \frac{1}{8} \frac{\beta^2}{a^2} \left( 1 + \frac{\beta^2}{a^2} \right) \begin{pmatrix} 3 & 1 \\ 1 & -1 \end{pmatrix} \quad (52)$$

Since we are interested in the implications of mimetic massive gravity for the growth and properties of large-scale structure in the late Universe, let us focus on subhorizon scales (i.e., $k^2 \gg a^2 H^2$) and assume the quasi-static (QS) approximation. In order to use this approximation, we first need to ensure that fluctuations in this régime are stable. Ignoring time variation in $a(t)$, which will be subdominant in the limit $k^2 \gg a^2 H^2$, and assuming solutions of the form $\vec{\chi} = \vec{\chi}_0 e^{i\omega t}$, the equations of motion following from the action (49) are

$$\left(-\omega^2 \mathbb{K} + \frac{k^2}{a^2} \mathbb{G} + m^2 \mathbb{M}\right) \vec{\chi} = 0. \quad (53)$$

We can then derive stability conditions from the dispersion relations, obtained by solving

$$\begin{align} 0 &= \det \left( -\omega^2 \mathbb{K} + \frac{k^2}{a^2} \mathbb{G} + m^2 \mathbb{M} \right) \nonumber \\ &= \frac{\omega^4}{a^2 + \beta^2} - \frac{\omega^2 k^2}{a^2 (3a^2 - \beta^2)} - \frac{5\omega^2 m^2 \beta^2}{8a^4} + \frac{k^2 m^2 \beta^2}{8a^6} \nonumber \\ &\quad + \frac{m^4 \beta^4 (a^2 + \beta^2)}{16a^8} \tag{54} \end{align}$$

for $\omega^2$.

The dispersion relations arising from eq. (54) are complicated, but simplify significantly in the limit $k \gg aH$ when we take into account the constraints (44) on $m/H_0$, which we obtained by requiring that the radiation and curvature densities generated by the mass term not exceed observational bounds. Consider replacing $m$ and $\beta$ in eq. (54) with the following two parameters,<fnref n="9" />

$$\epsilon_1 \equiv \left( \frac{m\beta}{k} \right)^2, \quad \epsilon_2 \equiv \left( \frac{m\beta^2}{ka} \right)^2. \quad (55)$$

We proceed to show that the bounds (44) imply that each of these is much smaller than unity on subhorizon scales for all times after matter-radiation equality.

<footnote n="9">⁹ To do this replacement, first replace $m \to \sqrt{\epsilon_1 \beta/k}$, and then replace any remaining factors of $\beta$ with $\beta \to \sqrt{\epsilon_2/\epsilon_1 a}$.</footnote>
For both $\epsilon_1$ and $\epsilon_2$ we can put upper bounds on the numerators and lower bounds on the denominators. Let us start with the numerators. For $\epsilon_1$, multiply each side of eq. (44) by $\beta$. We see there is a strict upper bound on the combination $m\beta$,

$$m\beta \le \sqrt{8\delta_K} H_0 \approx 0.15 H_0 \quad (56)$$

where we have taken $\delta_K \sim 0.003$ as a representative value. We can similarly find a bound on the numerator of $\epsilon_2$ by multiplying both sides of eq. (44) by $\beta^2$, finding

$$m\beta^2 \le 4\sqrt{\delta_r} H_0 \approx 10^{-2} H_0 \quad (57)$$

for $\delta_r \sim 10^{-5}$.

Now we move on to the denominators. The subhorizon limit is given by $k \gg aH$. For the sake of argument let us be conservative and assume that $k$ is only slightly subhorizon, $k/a \approx \mathcal{O}(1)H$.<fnref n="10" /> At any given time from matter-radiation equality to the present, where we can trust our analysis, the Hubble rate $H$ is related to its present-day value $H_0$ by $H = H_0\sqrt{\Omega_{\Lambda,0} + \Omega_{m,0}a^{-3}}$. Putting this together with the bounds we have derived on $m\beta$ and $m\beta^2$, we find

$$\epsilon_1 \lesssim \frac{0.02}{\Omega_{\Lambda,0}a^2 + \Omega_{m,0}a^{-1}} \ll 1, \quad (58)$$

$$\epsilon_2 \lesssim \frac{10^{-4}}{\Omega_{\Lambda,0}a^4 + \Omega_{m,0}a} \ll 1 \text{ for } z \lesssim 3000. \quad (59)$$

Note that while the upper bound on $\epsilon_1$ is always much smaller than unity for $0 < a \le 1$, the upper bound on $\epsilon_2$ in fact grows as $a^{-1}$ at early times. However, it grows slowly and has a factor of $10^{-4}$ to compete with, so that $\max(\epsilon_2)$ does not reach unity until $z \sim 3000$, right around matter-radiation equality. Therefore in principle there might be a handful of modes—right around the horizon scale and at the earliest moments of matter domination—for which terms going as $\epsilon_2$ affect the subhorizon dispersion relation, if $m\beta^2$ takes the largest value allowed by the constraints. We will continue to take $\epsilon_2 \ll 1$, with the understanding that if this particular situation is realized, then at those very early times we are only considering modes with $k \gtrsim 10aH$, for which $\epsilon_2$ is certainly smaller than unity.

Dropping terms subdominant in $\epsilon_1$ and $\epsilon_2$, the dispersion relation (54) becomes

$$0 \approx \frac{\omega^4}{a^2 + \beta^2} - \frac{\omega^2 k^2}{a^2(3a^2 - \beta^2)} + \frac{k^2 m^2 \beta^2}{8a^6}. \quad (60)$$

Solving for $\omega^2$, and again dropping terms subleading in $\epsilon_1 = (m\beta/k)^2$, we find the dispersion relations for our two modes,

$$\omega^2 \approx \frac{k^2}{a^2} \frac{a^2 + \beta^2}{3a^2 - \beta^2}, \quad (61)$$

$$\omega^2 \approx \frac{m^2 \beta^2}{8a^2} \left(3 - \frac{\beta^2}{a^2}\right). \quad (62)$$

Each of these implies the same stability condition,

$$\frac{\beta^2}{a^2} < 3. \quad (63)$$

This tells us that no matter what the value of $\beta$ is, our cosmological solutions are unstable at sufficiently high redshifts,

$$z > \sqrt{3}\beta^{-1} - 1. \quad (64)$$

This early time instability can however be safely pushed back to unobservably early times by taking the parameter $\beta$ to be sufficiently small.<fnref n="11" /> Because we are assuming matter and dark energy domination, we can trust our stability condition as far back as matter-radiation equality at $z_{eq} \approx 3400$. Demanding stability from $z_{eq}$ onward, we find a constraint on $\beta$,<fnref n="12" />

$$\beta \lesssim 5 \times 10^{-4}. \quad (65)$$

## 4.2. Cosmological tensor mass

Another possible cosmological bound on the parameters $m$ and $\beta$ comes from constraints on the graviton mass. The tightest bounds currently come from LIGO, $m_T \le 7.7 \times 10^{-23}$ eV [44].<fnref n="13" /> To compute the mass of tensor fluctuations on a cosmological background, we linearize the Einstein equation (3) around $g_{\mu\nu} = \tilde{g}_{\mu\nu} + h_{\mu\nu}$, with $\tilde{g}_{\mu\nu} = \text{diag}(-1, a^2\delta_{ij})$, $h_{00} = 0$, and $h_{ij}$ transverse and traceless, i.e., $h_{ii} = \partial_i h_{ij} = 0$. The Einstein equation is

$$\ddot{h}_{ij} + 3H\dot{h}_{ij} - \frac{\nabla^2}{a^2}h_{ij} + m_T^2 h_{ij} = 0 \quad (66)$$

with the tensor mass given by

$$m_T^2 = \frac{m^2}{2} \frac{\beta^2}{a^2} \left( 1 + \frac{\beta^2}{a^2} \right) \quad (67)$$

The structure of the Einstein equation is such that $m_T^2/m^2$ has to be a (quadratic) polynomial in $\beta^2/a^2$. What is non-trivial is that the degree-zero term in that polynomial cancels out, i.e., the expression for $m_T^2/m^2$ starts at order $\beta^2/a^2$. This means that gravitational waves propagating over cosmological distances (at low redshift, i.e., $a \sim \mathcal{O}(1)$) do not depend on $m$ alone; instead they involve the combinations $m\beta$ and $m\beta^2$ which, as we have seen, are strongly constrained by the cosmological background. In particular, recalling that $m^2\beta^2 \lesssim 10^{-2}H_0^2$ and $m^2\beta^4 \lesssim 10^{-4}H_0^4$, we see that $m_T$ at the present era is at most of order $10^{-1}H_0 \sim 10^{-34}$ eV, far below the LIGO bounds. Moreover, our stability condition (65) has no bearing on $m_T$. No matter how tiny $\beta$ is, the constraints (44) place a constant upper bound on $m\beta$, so that the smaller $\beta$ is, the larger $m$ is allowed to be, leaving $m_T \approx m\beta/(\sqrt{2}a)$ fixed. It is interesting to note that, without demanding that this model provides cosmic acceleration, the tensor mass is nevertheless forced to be smaller than the Hubble scale. Finally, we note that around a flat background, the tensor mass is simply $m$, so local tests of gravity might be able to place constraints on $m$ that are not possible with gravitational waves that propagate over cosmological distances.

## 4.3. Quasistatic limit

Finally, let us comment on the testability of mimetic massive gravity using near-future LSS surveys. We will find it convenient

<footnote n="10">¹⁰ Of course, the deeper in the subhorizon régime $k$ is, the smaller $\epsilon_1$ and $\epsilon_2$ become.</footnote>

<footnote n="11">¹¹ This is similar to massive bimetric gravity, which possesses an early-time instability that can be rendered safe in the limit where the ratio of the two Planck masses becomes small [43].</footnote>

<footnote n="12">¹² It is plausible that the result (63) holds, at least on an order-of-magnitude basis, through radiation domination as well (see, again, the example of biggravity [43]). In this case, we should demand that the instability be pushed back to before big bang nucleosynthesis, with $z_{BBN} \approx 3 \times 10^8$, which would imply a stronger constraint of $\beta \lesssim 10^{-8}$. We do not have much observational handle on the presumably radiation-dominated era before BBN, and therefore should not demand that the instability be absent then; indeed, a mild enough instability might have interesting consequences, such as the formation of primordial black holes.</footnote>

<footnote n="13">¹³ See Ref. [45] for a helpful summary of bounds on the graviton mass from a variety of experiments and observations.</footnote>
to work in Newtonian gauge, $B = E = 0$. Linearizing the Einstein equations (3), and leaving in a generic stress-energy tensor $T_{\mu\nu}$ for completeness, we obtain

$$6H^2\phi - \frac{2}{a^2}\partial^i\partial_i\psi + 6H\dot{\psi} = \frac{1}{M_{Pl}^2}\delta T_0^0 + 2\frac{\delta\lambda}{M_{Pl}^2} - \frac{m^2}{4}\frac{\beta^2}{a^2}\left(3 - \frac{\beta^2}{a^2}\right)\left(3\psi + \partial^i\partial_i\pi\right), \quad (68)$$

$$-2\partial_i(\dot{\psi} + H\phi) = \frac{1}{M_{Pl}^2}\delta T_i^0 + 2\frac{\bar{\lambda}}{M_{Pl}^2}\partial_i\pi^0 + \frac{m^2}{4}\left(3 - \frac{\beta^2}{a^2}\right)\left(\partial_i\pi^0 - \beta^2\partial_i\dot{\pi}\right), \quad (69)$$

$$6\left[\ddot{\psi} + 3H\dot{\psi} + H\Phi + (3H^2 + 2\dot{H})\phi\right] + \frac{2}{a^2}\partial^i\partial_i(\phi - \psi) = \frac{1}{M_{Pl}^2}\delta T_i^i - \frac{m^2}{4}\frac{\beta^2}{a^2}\left(1 + \frac{\beta^2}{a^2}\right)\left(3\psi + \partial^i\partial_i\pi\right), \quad (70)$$

$$\frac{1}{a^2}\partial^i\partial_j(\psi - \phi) = \frac{1}{M_{Pl}^2}\delta T_j^i + \frac{m^2}{2}\frac{\beta^2}{a^2}\left(1 + \frac{\beta^2}{a^2}\right)\partial^i\partial_j\pi, \quad (71)$$
$$i \neq j.$$

Moving to Fourier space, specializing to a pressureless fluid without anisotropic stress, and taking the quasistatic limit, $\dot{X} \sim H\dot{X} \sim H^2X \ll k^2X$ for any perturbation $X$, eqs. (68), (70) and (71) become

$$\frac{2k^2}{a^2}\psi = \frac{1}{M_{Pl}^2}(2\delta\lambda - \bar{\rho}\delta) - \frac{m^2}{4}\frac{\beta^2}{a^2}\left(3 - \frac{\beta^2}{a^2}\right)\left(3\psi - k^2\pi\right), \quad (72)$$

$$\frac{2k^2}{a^2}(\phi - \psi) = \frac{m^2}{4}\frac{\beta^2}{a^2}\left(1 + \frac{\beta^2}{a^2}\right)\left(3\psi - k^2\pi\right), \quad (73)$$

$$\frac{1}{a^2}(\phi - \psi) = -\frac{m^2}{2}\frac{\beta^2}{a^2}\left(1 + \frac{\beta^2}{a^2}\right)\pi, \quad (74)$$

where $\bar{\rho}$ and $\delta$ are the background density and overdensity of the dust component. Note that these are degenerate with the mimetic dark matter, as expected.

Combining these equations, we obtain the modified Poisson equation and the slip relation,

$$-k^2\psi = 4\pi G\mu(a, k)a^2(\delta\rho - 2\delta\lambda), \quad (75)$$

$$\psi = \eta(a, k)\phi, \quad (76)$$

where the modified-gravity parameters $\mu$ and $\eta$ are given by

$$\mu(a, k) = \frac{1}{1 + \frac{1}{2}\frac{m^2\beta^2}{k^2}\left(3 - \frac{\beta^2}{a^2}\right)}, \quad (77)$$

$$\eta(a, k) = \frac{1}{1 + \frac{1}{2}\frac{m^2\beta^2}{k^2}\left(1 + \frac{\beta^2}{a^2}\right)}. \quad (78)$$

These parametrize observable deviations from general relativity, in which $\mu = \eta = 1$.

The constraints we have already derived on $m$ and $\beta$ preclude $\mu$ and $\eta$ from deviating from unity at a level accessible to near-future observations. The stability constraint (63) requires the terms in parentheses to be $O(1)$, while the background constraint (56) sets $m^2\beta^2 \lesssim 0.02H_0^2$, so that

$$\mu - 1 \sim \eta - 1 \sim O\left(\frac{m^2\beta^2}{k^2}\right) \lesssim 10^{-2}\left(\frac{H_0}{k}\right)^2. \quad (79)$$

It is therefore highly unlikely that cosmological observations will be able to test this model against $\Lambda$CDM in the linear and sub-horizon régime.

## 5. Conclusions

In this Letter we have studied the first cosmological implications of the recently-proposed theory of mimetic massive gravity. We find that the theory is unable to self-accelerate without introducing a ghost. Its effects on Friedmann-Lemaître-Robertson-Walker cosmological backgrounds are to introduce effective radiation, curvature, and cosmological constant terms, as well as a dust-like mimetic dark matter component. We place constraints (44) on the theory parameters by demanding that the effective radiation and curvature terms be within observational bounds. In the ghost-free region of parameter space, $m^2 > 0$, the effective cosmological constant is negative-definite, so a separate dark energy sector, which we take to be a positive cosmological constant, is required to explain the late-time acceleration of the Universe.

We further studied the behavior of cosmological perturbations in the subhorizon, quasistatic limit. The model generically suffers from an instability at early times. However, since our analysis only included a pressureless dust component (in addition to a cosmological constant), the calculation can only be trusted as far back as matter-radiation equality. This allowed us to place a further constraint on the theory parameters by insisting that the instability be absent throughout matter domination. With these constraints, the deviations from $\Lambda$CDM in the linear, subhorizon régime are likely too small to be observable.

Not surprisingly, since this is a theory of massive gravity, it predicts massive tensor modes. We have calculated the tensor mass around cosmological backgrounds and found that, taking into account the constraints imposed by the cosmological background, this mass must be at least an order of magnitude below the Hubble scale, far outside the currently-available constraints on the graviton mass. Unlike other theories of massive gravity, in which the graviton mass is comparable to the Hubble scale in order to provide late-time acceleration, this bound on the graviton mass is solely due to the requirement that the effective radiation and curvature terms in the Friedmann equation not be too large.

What are the remaining prospects for cosmological tests of mimetic massive gravity? We emphasize that our analysis does not apply in two important régimes: horizon-size scales and nonlinear scales. One or both of these may possess signatures which could be used to distinguish mimetic massive gravity from $\Lambda$CDM, or otherwise to rule out additional regions of parameter space. One expects that nonlinear scales will require N-body simulations, while at horizon-size scales we cannot apply the quasistatic approximation and would need to solve the perturbation equations numerically, as in other theories of modified gravity [46]. For the latter, we note that the mass scales appearing in the action (49) for cosmological perturbations are not simply $m$, which can be arbitrarily large (in the limit of small $\beta$), but rather $m\beta$ and $m\beta^2$, which we have shown must both be at least an order of magnitude smaller than the Hubble scale. It therefore might be difficult for this theory to produce effects at horizon scales that are larger than cosmic variance. Note that scales $k \sim m\beta$ and $k \sim m\beta^2$ are super-horizon and therefore not observable.

## Acknowledgements

We are grateful to Marco Crisostomi, Sergei Dubovsky, Nick Kaiser, Riccardo Penco, Jeremy Sakstein, and Filippo Vernizzi for
helpful discussions. We thank the referee for useful comments on the draft. A.R.S. is supported by DOE HEP grants DOE DE-FG02-04ER41338 and FG02-06ER41449 and by the McWilliams Center for Cosmology, Carnegie Mellon University. V.V. is supported by a de Sitter PhD fellowship of the Netherlands Organization for Scientific Research (NWO). Y.A. is supported by LabEx ENS-ICFP: ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*. Y.A. also acknowledges support from the NWO and the Dutch Ministry of Education, Culture and Science (OCW), and also from the D-ITP consortium, a program of the NWO that is funded by the OCW.

## References

[1] A.H. Chamseddine, V. Mukhanov, J. High Energy Phys. 06 (2018) 060, arXiv: 1805.06283 [hep-th].

[2] A.H. Chamseddine, V. Mukhanov, J. High Energy Phys. 06 (2018) 062, arXiv: 1805.06598 [hep-th].

[3] A.H. Chamseddine, V. Mukhanov, J. High Energy Phys. 11 (2013) 135, arXiv: 1308.5410 [astro-ph.CO].

[4] M. Fierz, W. Pauli, Proc. R. Soc. Lond. A 173 (1939) 211.

[5] C. Isham, A. Salam, J. Strathdee, Phys. Rev. D 3 (1971) 867.

[6] D. Boulware, S. Deser, Phys. Rev. D 6 (1972) 3368.

[7] P. Creminelli, A. Nicolis, M. Papucci, E. Trincherini, J. High Energy Phys. 0509 (2005) 003, arXiv:hep-th/0505147 [hep-th].

[8] C. de Rham, G. Gabadadze, Phys. Rev. D 82 (2010) 044020, arXiv:1007.0443 [hep-th].

[9] C. de Rham, G. Gabadadze, A.J. Tolley, Phys. Rev. Lett. 106 (2011) 231101, arXiv: 1011.1232 [hep-th].

[10] S. Hassan, R.A. Rosen, Phys. Rev. Lett. 108 (2012) 041101, arXiv:1106.3344 [hep-th].

[11] S. Hassan, R.A. Rosen, A. Schmidt-May, J. High Energy Phys. 1202 (2012) 026, arXiv:1109.3230 [hep-th].

[12] S. Hassan, R.A. Rosen, J. High Energy Phys. 1204 (2012) 123, arXiv:1111.2070 [hep-th].

[13] S. Hassan, A. Schmidt-May, M. von Strauss, Phys. Lett. B 715 (2012) 335, arXiv: 1203.5283 [hep-th].

[14] S. Hassan, R.A. Rosen, J. High Energy Phys. 1202 (2012) 126, arXiv:1109.3515 [hep-th].

[15] K. Hinterbichler, R.A. Rosen, J. High Energy Phys. 07 (2012) 047, arXiv:1203. 5783 [hep-th].

[16] K. Hinterbichler, Rev. Mod. Phys. 84 (2012) 671, arXiv:1105.3735 [hep-th].

[17] C. de Rham, Living Rev. Relativ. 17 (2014) 7, arXiv:1401.4173 [hep-th].

[18] A.R. Solomon, Cosmology Beyond Einstein, Ph.D. thesis, Cambridge U., Cham, 2015, arXiv:1508.06859 [gr-qc].

[19] A. Schmidt-May, M. von Strauss, J. Phys. A 49 (2016) 183001, arXiv:1512.00021 [hep-th].

[20] E.A. Lim, I. Sawicki, A. Vikman, J. Cosmol. Astropart. Phys. 1005 (2010) 012, arXiv:1003.5751 [astro-ph.CO].

[21] L. Mirzagholi, A. Vikman, J. Cosmol. Astropart. Phys. 1506 (2015) 028, arXiv: 1412.7136 [gr-qc].

[22] H. van Dam, M.J.G. Veltman, Nucl. Phys. B 22 (1970) 397.

[23] V.I. Zakharov, JETP Lett. 12 (1970) 312; Pis'ma Zh. Eksp. Teor. Fiz. 12 (1970) 447.

[24] A. Higuchi, Nucl. Phys. B 282 (1987) 397.

[25] A. Vainshtein, Phys. Lett. B 39 (1972) 393.

[26] E. Babichev, C. Deffayet, Class. Quantum Gravity 30 (2013) 184001, arXiv:1304. 7240 [gr-qc].

[27] S.L. Dubovsky, J. High Energy Phys. 10 (2004) 076, arXiv:hep-th/0409124 [hep-th].

[28] D. Blas, D. Comelli, F. Nesti, L. Pilo, Phys. Rev. D 80 (2009) 044025, arXiv:0905. 1699 [hep-th].

[29] M. Lagos, M. Bañados, P.G. Ferreira, S. García-Sáenz, Phys. Rev. D 89 (2014) 024034, arXiv:1311.3828 [gr-qc].

[30] H. Motohashi, T. Suyama, K. Takahashi, Phys. Rev. D 94 (2016) 124021, arXiv: 1608.00071 [gr-qc].

[31] S. Dubovsky, T. Gregoire, A. Nicolis, R. Rattazzi, J. High Energy Phys. 03 (2006) 025, arXiv:hep-th/0512260 [hep-th].

[32] S. Endlich, A. Nicolis, R. Rattazzi, J. Wang, J. High Energy Phys. 04 (2011) 102, arXiv:1011.6396 [hep-th].

[33] G.R. Dvali, G. Gabadadze, M. Porrati, Phys. Lett. B 485 (2000) 208, arXiv:hep-th/0005016 [hep-th].

[34] C. Deffayet, Phys. Lett. B 502 (2001) 199, arXiv:hep-th/0010186 [hep-th].

[35] C. Deffayet, G.R. Dvali, G. Gabadadze, Phys. Rev. D 65 (2002) 044023, arXiv: astro-ph/0105068 [astro-ph].

[36] D. Gorbunov, K. Koyama, S. Sibiryakov, Phys. Rev. D 73 (2006) 044016, arXiv: hep-th/0512097 [hep-th].

[37] C. Charmousis, R. Gregory, N. Kaloper, A. Padilla, J. High Energy Phys. 10 (2006) 066, arXiv:hep-th/0604086 [hep-th].

[38] J. Khoury, J. Sakstein, A.R. Solomon, J. Cosmol. Astropart. Phys. 1808 (2018) 024, arXiv:1805.05937 [hep-th].

[39] N. Aghanim, et al., Planck, arXiv:1807.06209 [astro-ph.CO], 2018.

[40] X.-I. Chen, R.J. Scherrer, G. Steigman, Phys. Rev. D 63 (2001) 123504, arXiv: astro-ph/0011531 [astro-ph].

[41] R.H. Cyburt, B.D. Fields, K.A. Olive, T.-H. Yeh, Rev. Mod. Phys. 88 (2016) 015004, arXiv:1505.01076 [astro-ph.CO].

[42] D.J. Fixsen, Astrophys. J. 707 (2009) 916, arXiv:0911.1955 [astro-ph.CO].

[43] Y. Akrami, S.F. Hassan, F. Könnig, A. Schmidt-May, A.R. Solomon, Phys. Lett. B 748 (2015) 37, arXiv:1503.07521 [gr-qc].

[44] B.P. Abbott, et al., LIGO Scientific, VIRGO, Phys. Rev. Lett. 118 (2017) 221101, arXiv:1706.01812 [gr-qc]; Erratum: Phys. Rev. Lett. 121 (12) (2018) 129901.

[45] C. de Rham, J.T. Deskins, A.J. Tolley, S.-Y. Zhou, Rev. Mod. Phys. 89 (2017) 025004, arXiv:1606.08462 [astro-ph.CO].

[46] J. Enander, Y. Akrami, E. Mörtsell, M. Renneby, A.R. Solomon, Phys. Rev. D 91 (2015) 084046, arXiv:1501.02140 [astro-ph.CO].