# Massive and modified gravity as self-gravitating media

Guillermo Ballesteros, Comelli Denis, Pilo Luigi

► To cite this version:

Guillermo Ballesteros, Comelli Denis, Pilo Luigi. Massive and modified gravity as self-gravitating media. Physical Review D, 2016, 94, pp.124023. 10.1103/PhysRevD.94.124023. cea-01296404

HAL Id: cea-01296404

https://cea.hal.science/cea-01296404

Submitted on 1 Dec 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Massive and modified gravity as self-gravitating media

Guillermo Ballesteros,<fnref n="1" />,²* Denis Comelli,³,† and Luigi Pilo⁴,⁵,‡

¹Institut de Physique Théorique, Université Paris Saclay, CEA, CNRS, 91191 Gif-sur-Yvette, France

²Theory Division, CERN, 1211 Geneva, Switzerland

³Sezione di Ferrara, INFN, I-35131 Ferrara, Italy

⁴Dipartimento di Scienze Fisiche e Chimiche, Università di L'Aquila, I-67010 L'Aquila, Italy

⁵Laboratori Nazionali del Gran Sasso, INFN, I-67010 Assergi, Italy

(Received 11 June 2016; published 19 December 2016)

We study the effective field theory that describes the low-energy physics of self-gravitating media. The field content consists of four derivatively coupled scalar fields that can be identified with the internal comoving coordinates of the medium. Imposing SO(3) internal spatial invariance, the theory describes supersolids. Stronger symmetry requirements lead to superfluids, solids and perfect fluids, at lowest order in derivatives. In the unitary gauge, massive gravity emerges, being thus the result of a continuous medium propagating in spacetime. Our results can be used to explore systematically the effects and signatures of modifying gravity consistently at large distances. The dark sector is then described as a self-gravitating medium with dynamical and thermodynamic properties dictated by internal symmetries. These results indicate that the divide between dark energy and modified gravity, at large distance scales, is simply a gauge choice.

DOI: 10.1103/PhysRevD.94.124023

# I. INTRODUCTION AND OUTLINE

The cosmological data show that the Universe is under-
going a phase of accelerated expansion. See [1] for the
first conclusive evidence. Although a simple model of the
Universe based on a Friedmann-Lemaître-Robertson-
Walker (FLRW) metric in general relativity (GR) supple-
mented by a cosmological constant fits well all the
observations [2,3], the actual nature of the driving force
behind the acceleration remains still unclear. Upcoming
probes of the dynamics of the Universe [4] will be crucial to
shed light on this dark energy (DE). A good deal of
theoretical effort has been aimed in the last twenty years to
build consistent and compelling models of DE; see [5] for
reviews.<fnref n="1" /> What the vast majority of these efforts have in
common is the addition of some new degrees of freedom to
the dynamics of the Universe. On the one hand, these
hypothetical degrees of freedom have often been inter-
preted as extra components of the Universe beyond
baryonic and dark matter, photons and neutrinos. On the
other hand, they have also been commonly regarded as part
of gravity itself, modifying the behavior of GR at large

distances, in a way that is compatible with the current
acceleration of the Universe.

In this work we propose a symmetry-driven approach to
DE. By working with an action of four derivatively coupled
scalar fields, we show that DE can originate from a medium
with concrete mechanical and thermodynamic properties.
These four scalars are interpreted as comoving coordinates
of the medium and, moreover, they are the degrees of
freedom required to describe its low-energy physics. At the
core of this description of DE lie the symmetries of the
medium. They rule the mutual interactions of the four
scalsars and determine whether the medium behaves as a
perfect fluid, a superfluid, an elastic solid or a supersolid.

By using this large class of Lagrangians, we will argue
that the distinction that is often made between DE (as
some kind of matter) and modified gravity (as a true
modification of Einsten's gravity) is purely artificial.
Indeed both interpretations should be understood as
complementary points of view of a unique (hypothetical)
phenomenon. This fact becomes transparent by choosing
the adequate coordinates for each perspective, as we
explain below.

The four scalars can also be viewed as Stückelberg fields
that allow us to restore broken diffeomorphisms in models
of so-called massive gravity [7,8]. Following common
usage, we use the term massive gravity (MG) to refer to
any model of gravity for which nonderivative metric
perturbations appear squared. We recall, however, that
many such models do not contain a massive graviton (of
spin two). The introduction of these fields is the key for
interpreting massive (and general classes of modified)
gravity models as a self-gravitating media. Thanks to

<footnote>*guillermo.ballesteros@cea.fr</footnote>

<footnote>†comelli@fe.infn.it</footnote>

<footnote>‡luigi.pilo@aquila.infn.it</footnote>

<footnote n="1">¹However, it is worth stressing that no such model solves the problem of the smallness of the cosmological constant [6].</footnote>

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 3.0 License. Further distri-
bution of this work must maintain attribution to the author(s) and
the published article's title, journal citation, and DOI.</footnote>
diffeomorphism invariance, an adequate choice of the spacetime coordinates (called unitary gauge) renders trivial the dynamics of the Stückelberg fields, allowing us to draw a direct connection between the theory of self-gravitating media and MG.<fnref n="2" /> In particular, we describe in this work the matching to the general Hamiltonian construction of MG models that was studied in [11–13].

Most of such MG models are Lorentz breaking in the sense that, around Minkowski spacetime, the metric perturbations are just rotationally invariant in three dimensions. From the point of view of the self-gravitating medium interpretation of MG that we explore in this work, Lorentz breaking is simply triggered by the presence of the medium itself, which generically defines at least one preferred four-vector.

By construction, the general description of self-gravitating media in terms of four derivatively coupled scalars is a well-defined effective field theory (EFT). This allows a systematic exploration of DE in terms of symmetries, with a unified and suggestive interpretation of modified and MG models. This is the take-home message of this paper. The use of an EFT framework to describe modified gravity at large distances and the acceleration of the Universe, opens the possibility of confronting simultaneously a wide variety of models with cosmological observations in a simple way.<fnref n="3" />

The outline of the paper is the following. In Sec. II we start introducing the EFT of self-gravitating media, emphasizing the role of diffeomorphism invariance. In Sec. III we continue the presentation of the EFT by building the intuition that allows us to interpret the low-energy degrees of freedom as comoving coordinates. In Sec. IV we write the lowest-order scalars of the EFT according to the possible symmetries. In Sec. V we provide a discussion of the different kinds of media according to their symmetries and their field content. In Sec. VI the unitary gauge is employed to show how the EFT of self-gravitating media can be put in correspondence with a broad class of MG models. In Sec. VII we set the basis for a systematic study of the cosmology of these models by determining, for each kind of medium, the number of propagating degrees of freedom in a FLRW background. We present our conclusions in Sec. IX. Two appendixes are also included. In Appendix A we discuss briefly the conserved currents that appear in the theory. Finally, in Appendix B we give an alternative, shortcut derivation of the number of propagating degrees of freedom.

Through the paper we will use the signature convention $(-, +, +, +)$ and units such that $\hbar = c = 1$.

## II. DERIVATIVELY COUPLED SCALAR FIELDS AND DIFFEOMORPHISM INVARIANCE

Let us consider the action for a set of four scalar fields under diffeomorphisms $\Phi^A(x^\mu)$, $A = 0, 1, 2, 3$, whose mutual interactions obey shift symmetries,

$$ \Phi^A \rightarrow \Phi^A + c^A, \quad \partial_\mu c^A = 0. \qquad (2.1) $$

At lowest order in derivatives, the action for these fields depends on the kinetic blocks

$$ C^{AB} = g^{\mu\nu} \partial_{\mu} \Phi^{A} \partial_{\nu} \Phi^{B}, \qquad (2.2) $$

where $g_{\mu\nu}$ is the (four-dimensional) spacetime metric. Including the usual Einstein-Hilbert term the general form of the action at lowest order in derivatives is then given by

$$ S = \int d^4 x \sqrt{-g} [M_{pl}^2 R + U(C^{AB}) + L_m], \qquad (2.3) $$

where we define $M_{pl}^2 = 1/(16\pi G)$, $G$ being Newton's constant; $U$ is a smooth master function of the $4 \times 4$ matrix $C^{AB}$ and $L_m$ is the part of the Lagrangian describing matter fields (e.g. baryonic matter). The expression (2.3) is the leading contribution to the action of an EFT where shift symmetries enforce the fields $\Phi^A$ to be derivatively coupled.<fnref n="4" /> The scalars $C^{AB}$ are the only independent ones that can be built out of the four $\Phi^A$ at leading order (LO) in derivatives, i.e. with a single derivative acting on each $\Phi^A$. The next-to-leading order (NLO) in derivatives corresponds to operators that enter in the action with exactly two additional derivatives.<fnref n="5" /> These operators must then appear suppressed with respect to the LO ones by some energy scale $\Lambda$ squared. In principle, one would expect $\Lambda \sim U^{1/4}$, because this is the only scale of (2.3) other than $M_{pl}$ (and any masses coming from $L_m$).

At LO, we are also allowed to multiply the curvature scalar $R$ by an arbitrary function of $C^{AB}$. However, this function can always be eliminated through a conformal transformation of the metric,<fnref n="6" /> see [16]. After such a transformation, we obtain an action with the same structure of (2.3) plus some extra terms which can be disregarded at LO because they are of higher order in derivatives. Similarly, the Lagrangian $L_m$ can be multiplied by an arbitrary function of $C^{AB}$. In fact, denoting by $\psi_m$ the matter fields contained in $L_m$, according to the power

<footnote n="2">²An example of such media is the case of the perfect fluids that can be described with three Stückelberg fields, for which the invariance under internal volume preserving diffeomorphisms forbids a mass term for gravitons [SO(3) spatial tensor modes], see e.g. [9,10].</footnote>

<footnote n="3">³See also [14] for an approach with a similar spirit, where only a single new degree of freedom is considered and nonderivative couplings are allowed.</footnote>

<footnote n="4">⁴See [9] for a general discussion of this type of EFT with an arbitrary number of scalar fields.</footnote>

<footnote n="5">⁵See [15,16] for the NLO of a simpler case containing just three scalar fields.</footnote>

<footnote n="6">⁶We follow here the nomenclature of [17], but the reader must be aware that this is often called a Weyl transformation.</footnote>
counting this part of the action should actually be written as $f(C^{AB})L_m[h(C^{AB})g_{\mu\nu}, \psi_m]$, where $f$ and $h$ denote a set of arbitrary functions [16]. In this paper we will not study the effect of these couplings to matter fields, since our emphasis is on the interactions between the fields $\Phi^A$ and the metric $g_{\mu\nu}$. Therefore, in what follows we set $L_m$ to zero, postponing its study for future work.

Under these assumptions, the LO gravitational energy-momentum tensor (EMT) associated to the scalar fields $\Phi^A$, defined via the variation $\sqrt{-g}T_{\mu\nu}\delta g^{\mu\nu} = -2\delta \int d^4x \sqrt{-g}U$, is given by

$$T_{\mu\nu} = U g_{\mu\nu} - 2 \frac{\partial U}{\partial C^{AB}} \partial_{\mu} \Phi^A \partial_{\nu} \Phi^B, \quad (2.4)$$

where sums over the repeated indices A and B are implicit. Varying the action (2.3), the equations of motion (EOMs) of the scalars $\Phi^A$ read

$$\nabla_{\mu}(\nabla^{\mu}\Phi^{B}\partial U / \partial C^{AB}) = 0, \quad (2.5)$$

where $\nabla_{\mu}$ is the covariant derivative defined with the Levi-Civita connection of $g_{\mu\nu}$. Taking the covariant divergence of (2.4), and using that the Levi-Civita connection is torsionless, we obtain

$$\nabla^{\mu} T_{\mu\nu} = \nabla_{\nu} \Phi^{A} \nabla_{\mu} (\nabla^{\mu} \Phi^{B} \partial U / \partial C^{AB}). \quad (2.6)$$

Therefore, the covariant conservation of the EMT ($\nabla_{\mu}T^{\mu\nu} = 0$) can be inferred from the scalar's EOMs (2.5). Clearly, this result would also hold if $N$, the total number of scalars $\Phi^A$, were different from 4.

Generically, given some model of scalar fields coupled to gravity, the converse of the previous result is not true: the conservation of the EMT does not imply the EOMs of the scalars. However, in our case the two are equivalent under a specific condition. Let us consider (2.3) with $L_m = 0$ and assume (as we will do in most of this paper) that we have four scalars $\Phi^A$. If $\det(\partial_\mu \Phi^A) \neq 0$, it is possible to choose a system of coordinates such that

$$\partial_{\mu} \Phi^{A} = \delta_{\mu}^{A}. \quad (2.7)$$

Then, Eq. (2.6) becomes $\nabla^{\mu}T_{\mu\nu} = \delta_{\nu}^{A}\delta_{\mu}^{B}\nabla^{\mu}(\partial U / \partial C^{AB})$, which is nothing but (2.5) and shows that the EOMs of the scalars now follow from the conservation of the EMT. Since $\nabla^{\mu}T_{\mu\nu} = 0$ is ultimately due to diffeomorphism invariance, we can say that the equations of motion are a consequence of it, provided that $\det(\partial_{\mu}\Phi^A) \neq 0$.<fnref n="7" /> A more

direct way of proving it uses that (2.6) can be expressed as the product $Q \cdot q$, where $Q$ is a $4 \times 4$ matrix of components $\nabla_{\nu}\Phi^A$ and $q$ is a $4 \times 1$ matrix. Then, if $\det Q \neq 0$, the only solution of $Q \cdot q = 0$ is $q = 0$, which is (2.5).

It is well known that diffeomorphism invariance can be restored for any action of metric perturbations around some background. In four dimensions this can be done by introducing four scalars that are called Stückelberg fields.<fnref n="8" /> By construction, these scalar fields are coupled among themselves only derivatively. From this point of view, the action (2.3) with $L_m = 0$ is a good candidate for an EFT description of massive (and modified) gravity. The choice of $N = 4$ scalar fields $\Phi^A$ is then motivated by the number of spacetime dimensions we assume. As the Stückelberg fields acquire a nontrivial background configuration, such as e.g. (2.7), diffeomorphism invariance appears to be broken.

MG models (like the archetypical one of Fierz and Pauli [20]) are often presented in a way in which general covariance is explicitly broken (and in some cases Lorentz invariance too). In this work we exploit the fact that these models can generically be embedded in a covariant (and Lorentz invariant) action by the introduction of four scalar Stückelberg fields. As we have just explained, the apparent breaking of diffeomorphism invariance is just a consequence of choosing a background solution of the equations of motion for these fields. From this perspective, the action (2.3), together with a background satisfying (2.7), is an appropriate tool for the description of MG models, as we will discuss later in more detail.

The structure of (2.3) makes it convenient to assign dimensions of length to the fields $\Phi^A$. With this convention, the LO operators are dimensionless and the NLO ones have dimension of energy squared, which makes explicit the fact that they must appear suppressed by $\Lambda^2$. We can always express the Stückelberg fields $\Phi^A$ introducing new fields $\pi^A$ as follows:

$$\Phi^A(x^\mu) = x^A + \pi^A(x^\mu). \quad (2.8)$$

By construction, the fields $\pi^A$ inherit the power counting defined for the Stückelberg fields $\Phi^A$. Inserting (2.8) in the action (2.3), assuming that we can neglect NLO and higher-order operators, and imposing that $|\partial_{\mu}\pi^A| \ll 1$, we get an action for the $\pi^A$ that we can expand in powers of $\partial_{\mu}$ or, equivalently, in powers of $\pi^A$. With the condition $|\partial_{\mu}\pi^A| \ll 1$, the field redefinition (2.8) suggests that the

<footnote n="7">⁷An analogous argument can be applied if the number of scalar fields, $N$, is smaller than the number of spacetime dimensions, $D$. However, if $N > D$ there would not be enough freedom to eliminate all the scalars from the dynamics through (2.7) and so the EOMs cannot be equivalent the conservation of the EMT in that case.</footnote>

<footnote n="8">⁸The term Stückelberg field is commonly used for a field that allows to make explicit a (spontaneously broken) gauge symmetry. These fields and the associated mechanism (or “trick”), which we will use later on, are named after E.C.G. Stückelberg, who showed how a (complex) scalar field can make an Abelian vector field massive while preserving gauge invariance [18]. See [19] for a review with various applications.</footnote>
fields $\pi^A$ parametrize the deviation of the fields $\Phi^A$ from their background location $x^A$. Therefore, after the interpretation of $\Phi^A$ as coordinates in a medium, Eq. (2.8) tells us that $\pi^A$ may be understood as the degrees of freedom responsible for carrying sound waves in the medium itself, i.e. phonons which can be understood as the Goldstone bosons of broken translations, see e.g. [21,22]. We emphasize that these fields need not be small; it is only their derivatives that have to be under control for the validity of the derivative expansion. The fields $\pi^A$ are the four-dimensional generalization of the standard displacement field used in fluid dynamics, elasticity and strain theory. In the next section we start exploiting this connection, which we will phrase directly using the Stückelberg fields, using the language of the EFTs of solids, fluids and superfluids; different aspects of which have been studied in [9,10,15,16,23–33]. As a matter of fact, the basic idea of using three comoving coordinates to describe fluid dynamics is nothing but the earlier pull-back formalism, which has already been used to describe the dynamics of continuous media, see e.g. [34–39].

### III. SELF-GRAVITATING MEDIA

We split the scalars $\Phi^A$ into two distinct sets, $\{\Phi^0\}$ and $\{\Phi^1, \Phi^2, \Phi^3\}$, which are characterized as follows. The elements of $\{\Phi^1, \Phi^2, \Phi^3\}$ will be referred to as spatial Stückelberg fields, and we will use small latin letters $a, b, ...$ to index them. We require the $3 \times 3$ submatrix of (2.2) whose components are

$$B^{ab} \equiv C^{ab} \qquad (3.1)$$

to be positive definite, so that its three eigenvalues are positive. As we will soon see, the reason for this condition is allowing the inverse of $B^{ab}$ to define an Euclidean metric in $\mathbb{R}^3$. As anticipated in the previous sections, we can interpret the fields $\Phi^a$ as the (spatial) comoving coordinates of an internal space $\mathcal{F}$ that describes a continuum medium whose topology we assume to be $\mathbb{R}^3$. This picture—that has been known for a long time, see e.g. [35]—allows us to write an unconstrained variation principle for a fluid. A pedagogical review of this pull-back formalism is provided in [36] and various applications specific to cosmology can be found in [10,16,28,31,40–45]. Here, we will follow [10] to describe the basics of this construction.

We exploit the fact that, given a slicing and a threading of the spacetime (Lorentzian) manifold $\mathcal{M}$, we can formally write $\mathcal{M} = \mathcal{M}_T \times \mathcal{M}_S$, where $\mathcal{M}_T$ and $\mathcal{M}_S$ are respectively the sets of time, $\{x^0\}$, and space, $\{x^1, x^2, x^3\}$, coordinates on $\mathcal{M}$. By choice, $\mathcal{M}_S$ is endowed with a Riemannian metric. Provided that at any given time $x^0$, the condition $\det(\partial_i \Phi^a) \neq 0$ is satisfied, there exists an invertible map of $\mathcal{M}_S$ into the medium $\mathcal{F}$, which we assume to be smooth: $x^i \to \Phi^a(x^0, x^j)$. Conversely, $\mathcal{F}$ can be immersed

in $\mathcal{M}_S$ via the inverse map: $\Phi^a \to x^i(x^0, \Phi^a)$. While the first of these maps identifies the fluid (or in general, medium) element (labeled by $\Phi^a$) that sits at time $x^0$ at the point $x^i$; the second map describes the trajectory of the medium element $\Phi^a$ in spacetime [10]. This is nothing but the dual description of a continuum in terms of Lagrangian ($\Phi^a$) and Eulerian ($x^i$) coordinates.

The spacetime metric $g_{\mu\nu}$ induces a three-dimensional metric on $\mathcal{F}$, which is given by the matrix $B_{ab}$ defined as follows:

$$ds_{\mathcal{F}}^2 = g_{\mu\nu} \frac{\partial x^{\mu}}{\partial \Phi^a} \frac{\partial x^{\nu}}{\partial \Phi^b} d\Phi^a d\Phi^b \equiv B_{ab} d\Phi^a d\Phi^b. \qquad (3.2)$$

The induced metric $B_{ab}$ is the inverse matrix of $B^{ab}$ defined in (3.1), i.e. $B^{ab}B_{bc} = \delta_c^a$. The condition that $B^{ab}$ is positive definite ensures that the induced spatial metric on the medium $\mathcal{F}$ is Riemannian.

Since the fields $\Phi^a$ can be interpreted as comoving coordinates, they must remain unchanged along the fluid flow. This allows to define a (unique) four-velocity $u^\mu$, through the conditions

$$u^{\mu} \partial_{\mu} \Phi^{a} = 0, \qquad u^{\mu} u^{\nu} g_{\mu\nu} = -1, \qquad (3.3)$$

whose only solution is

$$u^{\mu} = - \frac{e^{\mu\nu a\beta}}{6b\sqrt{-g}} \epsilon_{abc} \partial_{\nu} \Phi^{a} \partial_{\alpha} \Phi^{b} \partial_{\beta} \Phi^{c}, \qquad (3.4)$$

where

$$b \equiv \sqrt{\det \mathbf{B}}, \qquad (3.5)$$

and $\mathbf{B}$ denotes the $3 \times 3$ matrix whose components are given by (3.1). In what follows, we will always use boldface capital letters for three-dimensional spatial matrices.

When pulled back into the spacetime manifold $\mathcal{M}$, the induced metric $B_{ab}$ becomes a projector $H_{\mu\nu}$ on the embedding of $\mathcal{F}$ in $\mathcal{M}$,

$$H_{\mu\nu} = B_{ab} \partial_{\mu} \Phi^{a} \partial_{\nu} \Phi^{b} \equiv g_{\mu\nu} + u_{\mu} u_{\nu}, \qquad H_{\mu\nu} u^{\nu} = 0. \qquad (3.6)$$

Only if $u_{[\alpha} \nabla_{\mu} u_{\nu]} = 0$ (where the brackets indicate full antisymmetrization of the three indices), Frobenius’s theorem guarantees that $u^\mu$ is hypersurface orthogonal. In that case there is a one-parameter family of hypersurfaces that are orthogonal (everywhere) to the four-velocity $u^\mu$, which thus defines a slicing of the spacetime manifold $\mathcal{M}$.

So far we have only dealt with the spatial Stückelberg fields, describing the structure that was used, e.g. in [10] to study (a certain class of) perfect fluids. We will now also deal with the temporal Stückelberg, $\Phi^0$, for which we impose the condition
$$X \equiv C^{00} < 0. \qquad (3.7)$$

This condition allows to define another four-velocity by means of the derivative $\partial_{\mu} \Phi^0$, which is thus constrained to be timelike,

$$V^{\mu} = -\frac{\partial^{\mu} \Phi^{0}}{\sqrt{-X}}, \qquad V^{\mu} V^{\nu} g_{\mu\nu} = -1. \qquad (3.8)$$

This defines a hypersurface at each point that is orthogonal to $V^{\mu}$ and has an induced metric

$$p_{\mu\nu} = g_{\mu\nu} + V_{\mu}V_{\nu}, \qquad p_{\mu\nu}V^{\nu} = 0. \qquad (3.9)$$

Since $V_{\mu}$ is the gradient of a scalar, its covariant vorticity $\omega_{\mu\nu} = p_{\mu}^{\alpha}p_{\nu}^{\beta}\nabla_{\alpha}V_{\beta}$ vanishes. In Minkowski and FLRW spacetimes, the four-vectors $u^{\mu}$ and $V^{\mu}$ coincide, up to a sign. However, if these spaces are slightly distorted, the two velocities generically differ from each other. Given the four scalar fields $\Phi^A$, the four-vectors $V^{\mu}$ and $u^{\mu}$ are the only independent four-vectors that are invariant under the group $SO(3)_s$ of internal spatial rotations and under the shift symmetries (2.1).

The condition (3.7) also allows for an interesting extension of the maps between $M_S$ and $\mathcal{F}$ that we mentioned earlier. In particular we can extend the medium space $\mathcal{F}$, spanned by the three spatial Stückelberg fields $\Phi^a$, to a space $\mathcal{F}_4$ spanned by the full set $\{\Phi^0, \Phi^a\}$. By construction, and in analogy to (3.2), the matrix $C_{AB}$, defined through

$$C_{AB} \equiv g_{\mu\nu} \frac{\partial x^{\mu}}{\partial \Phi^{A}} \frac{\partial x^{\nu}}{\partial \Phi^{B}}, \qquad (3.10)$$

is a Lorentzian metric in the extended medium space $\mathcal{F}_4$, and is the inverse of the $4 \times 4$ matrix of kinetic blocks introduced in (2.2). Hence, the map

$$M \to \mathcal{F}_4: x^{\mu} \to \Phi^A(x^{\mu}) \qquad (3.11)$$

carries the spacetime manifold onto the extended (self-gravitating) medium space, whereas the inverse

$$\mathcal{F}_4 \to M: \Phi^A \to x^{\mu}(\Phi^A) \qquad (3.12)$$

gives the location of each element of the extended medium (labeled by $\Phi^A$) on spacetime, as described by the coordinates $x^{\mu}$. In other words, the second of these two maps describes the world hypervolume of the extended medium $\mathcal{F}_4$ as it propagates on the spacetime $M$.

The role of the temporal Stückelberg $\Phi^0$ is therefore twofold. First, it allows us to define an independent timelike four-vector via (3.8). As we will later see, this will give us the possibility of defining two-constituent fluids, as required to describe superfluids (provided that the

appropriate symmetries are satisfied) [26]. And second, as we have just described, $\Phi^0$ can also be naturally interpreted as a comoving time coordinate in the medium space, which allows to conceive the self-gravitating media described by (2.3) as four-dimensional objects propagating in spacetime, extending the pull-back formalism to include the temporal dimension as well. Both pictures play important roles in the novel interpretation of MG and DE as self-gravitating media that we explore in this work.

In the next section we focus on the operators that can appear in the action (2.3) at LO in derivatives, and in Sec. V we will classify the different media that arise combining them.

## IV. SPATIAL SO(3) INVARIANCE AND ADDITIONAL SYMMETRIES

For simplicity, we will always impose a global internal spatial $SO(3)$ symmetry in the medium space, that we denote by $SO(3)_s$, so that all the operators in the action are invariant under

$$\Phi^a \to R^{ab} \Phi^b \quad \text{with} \quad R \in SO(3), \quad \partial_{\mu} R = 0. \qquad (4.1)$$

This basic symmetry will also allow us to make direct contact with the general models of massive gravity classified in [13], as we will discuss in Sec. VI.

To construct a LO complete set of independent $SO(3)_s$-invariant scalars we start with the traces

$$\tau_n = \mathrm{Tr}(\mathbf{B}^n), \qquad n = 1, 2, 3. \qquad (4.2)$$

We recall that the $SO(3)_s$-invariant determinant $b = \sqrt{\det \mathbf{B}}$, defined already in (3.5), is not independent of these traces because it can be written as $b^2 = (\tau_1^3 - 3\tau_1\tau_2 + 2\tau_3)/6$.

Introducing the $3 \times 3$ matrix $\mathbf{Z}$ of components

$$Z^{ab} = C^{0a}C^{0b}, \qquad (4.3)$$

we can define the scalars

$$y_n = \mathrm{Tr}(\mathbf{B}^n \cdot \mathbf{Z}), \qquad n = 0, 1, 2, 3, \qquad (4.4)$$

which are also invariant under $SO(3)_s$. The traces of $\mathbf{Z}^2$ and $\mathbf{Z}^3$ (and higher powers) do not need to be considered once $y_0$ is included because they can be written as powers of the latter.

In addition, since the four-vectors $\partial_{\mu} \Phi^0$ and $u_{\mu}$ are invariant under $SO(3)_s$, there are two other scalars at LO that are invariant under this group. One of them is $X = C^{00}$, which we already introduced in (3.7), and the other is

$$Y \equiv u^{\mu} \partial_{\mu} \Phi^{0}. \qquad (4.5)$$
This closes the set of independent LO $SO(3)_s$-invariant operators, which is therefore composed of nine elements and reads

$$ \mathcal{O}_{\text{LO}} = \{X, Y, \tau_1, \tau_2, \tau_3, y_0, y_1, y_2, y_3\}. \quad (4.6) $$

We point out that the determinant of the 4 × 4 matrix of components $C^{AB}$, defined in (2.2), is also invariant under $SO(3)_s$, but it is not independent of this set because it can be expressed as

$$ \sqrt{|\det(C^{AB})|} = Yb. \quad (4.7) $$

What is interesting of this operator is that is closely related to the pull back of the medium volume element in spacetime. Indeed

$$ d^4\Phi = \sqrt{-gb}Yd^4x. \quad (4.8) $$

This expression allows to pass from Eulerian to Lagrangian coordinates (and vice versa) in four dimensions, extending the three-dimensional treatment of these coordinates that was discussed in [10].

## A. Additional symmetries at LO

We can further constrain the form of the action (2.3) by requiring additional symmetries, leading to different types of self-gravitating media. Two examples are elastic solids and perfect fluids, which are derived from symmetries imposed on the spatial sector of the Stückelberg fields. We will see in Sec. V that by imposing symmetries that mix the spatial and temporal Stückelberg fields, new kinds of materials arise. For reference, we summarize in Table I the notation for the operators that we use in this paper. Let us now make a list of the additional symmetries of interest for us,

(i) Imposing volume-preserving diffeomorphism invariance of the spatial sector,

$$ V_s \text{Diff}: \Phi^a \rightarrow \Psi^a(\Phi^b), \quad \det\left(\frac{\partial\Psi^a}{\partial\Phi^b}\right) = 1, \quad (4.9) $$

the scalars $\tau_n$ have to combine in the action into the determinant operator $b$ and, in addition, the scalars $y_n$ are automatically forbidden. Therefore $b, X$ and $Y$ are the only operators allowed at LO. This symmetry was identified in [9] in the context of the EFT formulation of perfect fluids containing only $b$.

We will also consider a series of symmetries that were proposed in [7].

(i) If we impose that the action has to be invariant under the transformations

$$ \Phi^a \rightarrow \Phi^a + f^a(\Phi^0), \quad (4.10) $$

<table><caption>TABLE I. Summary of the various LO scalar operators appearing in this paper. Unless it is explicitly indicated otherwise, we use greek letters for spacetime indices, capital latin letters (A, B, ..., = 0, 1, 2, 3) for indices in the internal spacetime of the medium and small latin letters (a, b, ... = 1, 2, 3) for spatial indices of the medium. Boldface latin capital letters, such as **B**, represent 3 × 3 matrices.</caption><thead><tr><th>Operator</th><th>Definition</th></tr></thead><tbody><tr><td>$C^{AB}$</td><td>$g^{\mu\nu} \partial_{\mu} \Phi^A \partial_{\nu} \Phi^B$, A, B = 0, 1, 2, 3</td></tr><tr><td>$B^{ab}$</td><td>$g^{\mu\nu} \partial_{\mu} \Phi^a \partial_{\nu} \Phi^b$, a, b = 1, 2, 3</td></tr><tr><td>$Z^{ab}$</td><td>$C^{a0} C^{b0}$</td></tr><tr><td>$X$</td><td>$C^{00}$</td></tr><tr><td>$W^{ab}$</td><td>$B^{ab} - Z^{ab} / X$</td></tr><tr><td>$b$</td><td>$\sqrt{\det \mathbf{B}}$</td></tr><tr><td>$Y$</td><td>$u^{\mu} \partial_{\mu} \Phi^0$</td></tr><tr><td>$y_n$</td><td>$\text{Tr}(\mathbf{B}^n \cdot \mathbf{Z})$, n = 0, 1, 2, 3</td></tr><tr><td>$\tau_n$</td><td>$\text{Tr}(\mathbf{B}^n)$, n = 1, 2, 3</td></tr><tr><td>$w_n$</td><td>$\text{Tr}(\mathbf{W}^n)$, n = 0, 1, 2, 3</td></tr><tr><td>$\mathcal{O}_{\alpha\beta n}$</td><td>$(X/Y^2)^\alpha (y_n/Y^2)^\beta$, $\alpha, \beta \in \mathbb{R}$</td></tr><tr><td>$\mathcal{O}_\alpha$</td><td>$(X/Y^2)^\alpha$, $\alpha \in \mathbb{R}$</td></tr></tbody></table>

for arbitrary $f^a$, we find that at LO the action will generically depend only on $X$ and on the matrix of components

$$ W^{ab} \equiv B^{ab} - \frac{Z^{ab}}{X}. \quad (4.11) $$

Then, the invariance under $SO(3)_s$ further enforces $W^{ab}$ to appear in the action only through the traces

$$ w_n \equiv \text{Tr}(\mathbf{W}^n), \quad (4.12) $$

which are nonlinear combinations of the scalars $X$, $\tau_n$ and $y_n$,

$$ \begin{aligned} w_1 &= \tau_1 - \frac{y_0}{X}, & w_2 &= \tau_2 - 2\frac{y_1}{X} + \frac{y_0^2}{X^2}, \\ w_3 &= \tau_3 - 3\frac{y_2}{X} + 3\frac{y_0 y_1}{X^2} - \frac{y_0^3}{X^3}. & & \end{aligned} \quad (4.13) $$

Under the transformation (4.10), $u^\mu \partial_\mu \Phi^a \rightarrow u^\mu \partial_\mu \Phi^a + Y df^a / d\Phi^0$, because $u^\mu$ is invariant. Therefore, the operator $Y$ measures the failure of $\Phi^a$ to remain comoving along $u^\mu$ when such a transformation is applied. Besides, as in the case of $V_s \text{Diff}$, at NLO the same operators with derivatives of $u^\mu$ are also allowed under (4.10).

(ii) If the action is invariant under

$$ \Phi^0 \rightarrow \Phi^0 + f(\Phi^a), \quad (4.14) $$

the only LO operators that are allowed are $Y$ and the traces $\tau_n$. Clearly, $u^\mu$ is invariant under this
<table><caption>TABLE II. Summary of local symmetries in material spacetime and the corresponding invariant scalar operators. Invariance under $SO(3)_s$ and shift symmetries are assumed by default in all cases.</caption><thead><tr><th colspan="3">Four-dimensional media</th></tr><tr><th>Symmetries of the action</th><th>LO scalar operators</th><th>Type of medium</th></tr></thead><tbody><tr><td>$SO(3)_s$ & $\Phi^A \to \Phi^A + f^A$, $\partial_\mu f^A = 0$</td><td>$X, Y, \tau_n, y_n$</td><td>Supersolids</td></tr><tr><td>$\Phi^a \to \Phi^a + f^a(\Phi^0)$</td><td>$X, w_n$</td><td></td></tr><tr><td>$\Phi^0 \to \Phi^0 + f(\Phi^a)$</td><td>$Y, \tau_n$</td><td></td></tr><tr><td>$\Phi^0 \to \Phi^0 + f(\Phi^0)$</td><td>$\tau_n, w_n, O_{\alpha\beta n}$</td><td></td></tr><tr><td>$\Phi^a \to \Phi^a + f^a(\Phi^0)$ & $\Phi^0 \to \Phi^0 + f(\Phi^0)$</td><td>$w_n$</td><td></td></tr><tr><td>$V_s\text{Diff}: \Phi^a \to \Psi^a(\Phi^b)$, $\det|\partial\Psi^a/\partial\Phi^b| = 1$</td><td>$b, Y, X$</td><td>Superfluids</td></tr><tr><td>$\Phi^0 \to \Phi^0 + f(\Phi^0)$ & $V_s\text{Diff}$</td><td>$b, O_\alpha$</td><td></td></tr><tr><td>$\Phi^0 \to \Phi^0 + f(\Phi^a)$ & $V_s\text{Diff}$</td><td>$b, Y$</td><td>Perfect fluid</td></tr><tr><td>$\Phi^A \to \Psi^A(\Phi^B)$, $\det|\partial\Psi^A/\partial\Phi^B| = 1$</td><td>$bY$</td><td>Perfect fluid with $\rho + p = 0$</td></tr></tbody></table>

<table><caption>TABLE III. Summary of self-gravitating media with reduced dimensionality, i.e. with less than four Stückelberg fields.</caption><thead><tr><th colspan="3">Media with reduced internal dimensionality</th></tr><tr><th>Symmetries of the action</th><th>LO scalar operators</th><th>Type of medium</th></tr></thead><tbody><tr><td>$SO(3)_s$ & $\Phi^a \to \Phi^a + c^a$, $\partial_\mu c^a = 0$</td><td>$\tau_n$</td><td>Solid</td></tr><tr><td>$V_s\text{Diff}$</td><td>$b$</td><td>Perfect fluids</td></tr><tr><td>$\Phi^0 \to \Phi^0 + c^0$, $\partial_\mu c^0 = 0$</td><td>$X$</td><td></td></tr></tbody></table>

transformation and therefore all the NLO operators mentioned above constructed from derivatives of $u^\mu$ also respect this transformation.

If we combine the symmetries (4.14) and (4.9) the master function $U$ can only depend on $Y$ and $b$ and the resulting medium has the EMT of a perfect fluid, as we will later see.

(iii) Finally, if the action is invariant under

$$ \Phi^0 \to \Phi^0 + f(\Phi^0), \qquad (4.15) $$

it is constrained at LO to be a function of the traces $\tau_n, w_n$ and the scalars

$$ O_{\alpha\beta n} = \left(\frac{X}{Y^2}\right)^\alpha \left(\frac{y_n}{Y^2}\right)^\beta. \qquad (4.16) $$

The operators $O_{\alpha\beta n}$ present some peculiarities with respect to the ones that we have encountered so far. From the point of view of the symmetries, $\alpha, \beta \in \mathbb{R}$, which implies that there must be an infinite (and uncountable) amount of these operators. However, since all $y_n$ vanish at zero order for diagonal metric backgrounds, a well-defined perturbative expansion seems to require $\beta$ to be a positive integer. Even then, there is still an infinite number of these operators at LO, which means that the master function $U$ cannot be arbitrary if the action is to be finite. The same

difficulty is also manifest on the associated Noether currents, see Appendix A 2.

If in addition to (4.15), we impose also (4.9), the operators (4.16) get restricted to the subset with $\beta = 0$. For later convenience we will denote them as follows:

$$ O_\alpha = \left(\frac{X}{Y^2}\right)^\alpha. \qquad (4.17) $$

If instead we combine (4.15) and (4.10), the operators $w_n$ get selected.

In Table II we summarize the different LO scalar operator content according to the internal symmetries of the action and the physical interpretation of the resulting system, which we discuss in the next section. In addition, media with reduced internal dimensionality are possible as well; see Table III. For instance, if the spatial Stückelberg fields $\Phi^a$ are not present in the LO action, the only operator at that order would be $X$.

## V. LAGRANGIANS FOR SELF-GRAVITATING MEDIA AT LO

Given the field content, the metric $g_{\mu\nu}$ and the scalar fields $\Phi^A$, it is possible to study in a systematic way the models resulting from imposing certain symmetries in addition to $SO(3)_s$ and internal shifts. Each class of
symmetries corresponds to a medium with specific mechanical and thermodynamic properties, which can be obtained from the EMT. For all the media we consider we require spacetime diffeomorphism invariance. Then, the most general action at LO for a medium described by the four Stückelberg fields $\Phi^A$ can be constructed in terms of the scalar invariants $X, Y, \tau_n$ and $y_n$,

$$S = M_{pl} \int d^4 x \sqrt{-g} R + \int d^4 x \sqrt{-g} U(X, Y, \tau_n, y_n) \quad (5.1)$$

and the corresponding EMT is given by

$$T_{\mu\nu} = U g_{\mu\nu} - 2 \frac{\partial U}{\partial C^{AB}} \partial_{\mu} \Phi^A \partial_{\nu} \Phi^B = U g_{\mu\nu} - 2 \sum_k U_{O_k} \frac{\partial O_k}{\partial g^{\mu\nu}}, \quad (5.2)$$

where $O_k$ are the nine scalar LO operators appearing in (5.1) and we use the notation $U_{O_k} = \partial U / \partial O_k$. Their partial derivatives with respect to the inverse spacetime metric are

$$\begin{aligned}
\frac{\partial Y}{\partial g^{\mu\nu}} &= -\frac{Y}{2} u_{\mu} u_{\nu}, & \frac{\partial X}{\partial g^{\mu\nu}} &= -X \mathcal{V}_{\mu} \mathcal{V}_{\nu}, \\
\frac{\partial \tau_n}{\partial g^{\mu\nu}} &= n \partial_{\mu} \Phi \cdot \mathbf{B}^{n-1} \cdot \partial_{\nu} \Phi, & \frac{\partial b}{\partial g^{\mu\nu}} &= \frac{b}{2} H_{\mu\nu}, \\
\frac{\partial y_n}{\partial g^{\mu\nu}} &= \sum_{m=1}^{n} \partial_{\mu} \Phi \cdot \mathbf{B}^{n-m} \cdot \mathbf{Z} \cdot \mathbf{B}^{m-1} \\
& \qquad \cdot \partial_{\nu} \Phi - 2 \sqrt{-X} (C \cdot \mathbf{B}^n \cdot \partial_{(\mu} \Phi) \mathcal{V}_{\nu}), & (5.3)
\end{aligned}$$

where the indices $\mu$ and $\nu$ are symmetrized in the second term of the last expression. We have included in (5.3) the operator $b$, that can be written as a combination of the three $\tau_n$, for later convenience. The dot ($\cdot$) represents the standard three-dimensional matrix product and we have introduced the notation

$$C^i = C^{0i}. \quad (5.4)$$

In the expressions (5.3), $C$ and $\partial_{\mu} \Phi$ have to be understood as a $1 \times 3$ matrices of components $C^i$ and $\partial_{\mu} \Phi^i$, respectively. It is also useful to note that

$$\frac{\partial u^{\alpha}}{\partial g^{\mu\nu}} = -\frac{u^{\alpha}}{2} u_{\mu} u_{\nu}, \qquad \frac{\partial \mathcal{V}^{\mu}}{\partial g^{\alpha\beta}} = -\frac{\mathcal{V}^{\mu}}{2} \mathcal{V}_{\alpha} \mathcal{V}_{\beta}. \quad (5.5)$$

In what follows we list the media that arise according to the symmetries that we have discussed and, also, according to the assumptions on the low-energy field content. We focus first on the systems that contain only a partial subset of the Stückelberg fields, which we call temporal or spatial media. These are the solids $U(\tau_n)$, and two different types of perfect fluids, $U(b)$ and $U(X)$, the first of which is a special case of the solids. Then we will move on to the

media that at low energy require all the Stückelberg fields. As discussed below, the properties of the most generic class of them have been argued to define (nonrelativistic) supersolids [46]. A subclass of these (with an enhanced symmetry group: internal three-volume preserving diffeomorphisms $V_sDiff$) would then describe superfluids. Other media, defined by specific additional symmetries, will also be mentioned, see Table II. The connection with modified and, especially, massive gravity models will be studied in Sec. VI. Then, in Sec. VII we will discuss the number of propagating degrees of freedom in FLRW, which strongly depends on the symmetries beyond (2.1) and (4.1) that are imposed.

### A. Media with reduced internal dimensionality

We start by describing three special cases for which some of the four Stückelberg fields are missing from the LO action, while the shifts symmetries [and eventually while $SO(3)_s$] are respected. We refer to these media as temporal or spatial depending on whether the missing fields are the spatial or temporal Stückelbergs, respectively. Specifically, we will discuss two kinds of spatial media (solids and perfect fluids) and the only kind of temporal medium (an irrotational perfect fluid) that exists at LO with our symmetry assumptions. A peculiarity of the spatial and temporal media is that there is no symmetry that forbids one subset of Stückelberg fields at LO, but allows them at a higher order in derivatives. This means that the defining restriction on the low-energy field content of these models must be taken as an assumption, which in principle cannot be justified from the point of view of the symmetries.

#### 1. $U(\tau_1, \tau_2, \tau_3)$: Solid

Let us first consider the possibility that the only operators relevant at LO are the traces $\tau_n$. Then, the action contains only the three spatial Stückelberg fields $\Phi^a$. This necessarily assumes that no other fields are needed to describe the system at sufficiently low energies. As we just mentioned, there is no symmetry that forbids $\Phi^0$ at LO but allows it at some higher order in derivatives. Notice that by looking at Table II, it would seem that combining $\Phi^0 \to \Phi^0 + f(\Phi^0)$ and $\Phi^0 \to \Phi^0 + f(\Phi^a)$, the operators $\tau_n$ are selected. However, it is clear that this cannot be considered a true symmetry of an action based on $U(\tau_n)$ since $\Phi^0$ is not part of $\tau_n$.

In addition, given that the general action (5.1) lacks potential terms and it generically mixes all fields in a democratic way, it is difficult to imagine how $\Phi^0$ could be (on average) frozen at some constant value only at LO, if it is not neglected also at higher orders. Therefore, it seems that if $\Phi^0$ does not appear at LO, it should be completely absent from the action and this must be imposed as an assumption on the physical nature of the system.
The solid is thus the most general medium that we can construct at LO with only the three spatial Stückelberg fields, respecting the basic symmetries that we assumed (SO(3)$_s$ and shift invariance). The corresponding EMT can be written as

$$ T_{\mu\nu} = U g_{\mu\nu} - 2U_{ab}\partial_{\mu}\Phi^{a}\partial_{\nu}\Phi^{b}, \\ \text{where } U_{ab} = U_{\tau_1}\delta^{ab} + 2U_{\tau_2}B^{ab} + 3U_{\tau_3}B^{ac}B^{cb}. \quad (5.6) $$

Such a medium can be interpreted as the relativistic generalization of an elastic material, in the sense that this kind of media can support compressing pressures; see e.g. [47]. In standard elasticity theory in a flat space spacetime [48], the stress state is described by the spatial internal inverse metric $B^{ab}$, so that the medium is relaxed (not stressed) if $B^{ab} = \delta^{ab}$.

The object $B^{ab}$ is the covariant generalization of the standard three-dimensional left Cauchy-Green deformation tensor, and thus $\partial\Phi^a/\partial x^i$ is the deformation gradient tensor, see e.g. [49]. This type of medium, which has been commonly called solid (see e.g. [21]) or elastic solid (to distinguish it from a perfect fluid) has been applied in astrophysics and cosmology for the description of the dynamics of the interior of neutron stars [50], in a proposal for a dark matter model [51] and also for constructing a model of primordial inflation [28,40].

Solids posses two features, due to the independent operators $\tau_n$, that make them different from perfect fluids—which we discuss below—and interesting for certain applications in the context of cosmology. First, they exhibits anisotropic stress, due to the pressure perturbations being in general dependent on the spatial direction. And second, the spatial trace operators $\tau_n$ are the only ones among the nine LO scalars operators (4.6) that generate a mass for the (spin two) graviton, as we will see in Sec VI. Clearly, these two properties are not exclusive of the solids with EMT (5.6)—see Table II—but these are the simplest systems that display them.

Although it is not immediate to identify a four-velocity from the EMT (5.6), the natural choice is $u^\mu$, defined in (3.4), since this is the only SO(3)$_s$-invariant four-vector that is available when only the three spatial Stückelberg fields are present. With this choice of frame, the energy density and the pressure of the solid, defined as usual through the projections $\rho = T^{\mu\nu}u_\mu u_\nu$ and $3p = T^{\mu\nu}H_{\mu\nu} = T^{\mu\nu}(g_{\mu\nu} + u_\mu u_\nu)$, are given by

$$ \rho = -U, \quad p = U - \frac{2}{3} \sum_{n=1}^{3} n U_{\tau_n} \tau_n. \quad (5.7) $$

The four-velocity $u^\mu$ corresponds to the energy frame—sometimes called rest frame—of the solid, which for an arbitrary EMT is defined (if it exists) as the four-vector $U_{(r)}^\mu$ that solves $(g^{\mu\nu} + U_{(r)}^\mu U_{(r)}^\nu)U_{(r)}^\gamma T_{\nu\gamma} = 0$, see e.g. [10]. Then, the EMT (5.6) can be written as follows:

$$ T_{\mu\nu} = pg_{\mu\nu} + (\rho + p)u_{\mu}u_{\nu} + \pi_{\mu\nu}, \quad (5.8) $$

where the nonzero anisotropic stress $\pi_{\mu\nu}$ can be obtained, for instance, from the difference of (5.8) and (5.6). The expression (5.8) is what in cosmology is usually referred to as an imperfect fluid.

The symmetric matrix $B^{ab}$ can be decomposed in terms of its eigenvectors $\zeta_{(n)}^a$ and (real) eigenvalues $\lambda_n$, $n = 1, 2, 3$ as $B^{ab} = \sum_n \lambda_n \zeta_{(n)}^a \zeta_{(n)}^b$. Then, we can define the projectors $P_{(n)}^a = \zeta_{(n)}^a \zeta_{(n)}^b$ satisfying $P_{(n)} \cdot P_{(m)} = \delta_{mn}P_{(n)}$. The eigenvectors are mutually orthogonal with respect to the three-dimensional metric $\delta_{ab}$, i.e. $\delta_{mn} = \zeta_{(m)}^a \zeta_{(n)}^b \delta_{ab}$. By inspection, $U_{ab}$ is also diagonalized by $\zeta_{(m)}^a$ as follows:

$$ U_{ab} = \delta_{ac}\delta_{bd} \sum_{n=1}^{3} \tilde{\lambda}_{n} P_{(n)}^{cd}, \quad (5.9) $$

where $\tilde{\lambda}_n = U_{\tau_1} + 2U_{\tau_2}\lambda_n + 3U_{\tau_3}\lambda_n^2$. Therefore, the EMT (5.6) can then be cast in the following form:

$$ \begin{aligned} T_{\mu\nu} &= U g_{\mu\nu} - 2 \sum_{n=1}^{3} \tilde{\lambda}_n \lambda_n \zeta_{\mu}^{(n)} \zeta_{\nu}^{(n)}, \\ \zeta_{\mu}^{(n)} &= \lambda_n^{-1/2} \partial_{\mu} \Phi^a \zeta_{(n)}^b \delta_{ab}. \end{aligned} \quad (5.10) $$

Notice that $u^\mu \zeta_{\mu}^{(n)} = 0$ and $g^{\mu\nu} \zeta_{\mu}^{(m)} \zeta_{\nu}^{(n)} = \delta^{mn}$ and thus $\{u_\mu, \zeta_{\nu}^{(m)}\}$ form an orthonormal tetrad. The elastic solid has principal pressures given by

$$ p_n \equiv T_{\alpha\beta} \zeta^{(n)}_{\alpha} \zeta^{(n)}_{\beta} = U - 2 \tilde{\lambda}_n \lambda_n, \quad n = 1, 2, 3. \quad (5.11) $$

The thermodynamics of solids can be studied following analogous lines to those given below and in [52] for the special case of perfect fluids.

## 2. U(b): Perfect fluid

If the SO(3)$_s$ symmetry of the solid is enlarged to the invariance under volume preserving spatial diffeomorphisms $V_s$Diff, see (4.9), the trace operators $\tau_n$ must combine into the determinant $b^2$ defined in (3.5), constraining further the properties of the medium. Indeed, we recall that $b^2$ can be expressed as a function of $\tau_n$ using $b^2 = (\tau_1^3 - 3\tau_1\tau_2 + 2\tau_3)/6$. Then, under the assumption that only the spatial Stückelberg fields are present and imposing $V_s$Diff, the master function $U$ that determines the Lagrangian of the system depends only on $b$ at LO. However, if the restriction on the field content is relaxed, allowing also for the presence of the temporal Stückelberg $\Phi^0$, the operators $Y$ and $X$ should be included as well under at lowest order in derivatives, see Table II. In this section we will focus on the special case $U(b)$ and we will later move
progressively to the general structure of the media that respect $V_s$Diff.

In the case $U(b)$ the EMT reads

$$T_{\mu\nu} = pg_{\mu\nu} + (\rho + p)u_{\mu}u_{\nu}, \quad \rho = -U, \quad p = U - bU_b, \qquad (5.12)$$

with $u^\mu$ being defined in (3.4). In comparison to the general solid (5.8), the symmetry $V_s$Diff removes the anisotropic stress $\pi_{\mu\nu}$, thus restricting the medium to be a perfect fluid. In addition, it also sets to zero the (spin two) graviton mass.

The application of this type of systems for describing cosmological perfect fluids (and specifically their perturbations) was studied in [10]. For instance, if the fluid has a barotropic equation of state of the form $p = w\rho$, with constant $w$, the master function is of the form

$$U \propto b^{1+w}. \qquad (5.13)$$

Hence, the background dynamics of the $\Lambda$CDM universe can be modeled with

$$U \propto (\Omega_r(b^{4/3} - 1) + \Omega_m(b - 1) + 1), \qquad (5.14)$$

with $\Omega_r \sim 10^{-4}$ and $\Omega_m \sim 0.26$, to satisfy the current constraints on the radiation and (cold dark) matter densities [3].

In general, the perturbations of a $U(b)$ perfect fluid can be studied decomposing the spatial phonon fields $\pi^i$, $i = 1, 2, 3$—introduced in (2.8) for the general four-dimensional case—into longitudinal and transverse polarizations,

$$\pi^i = \pi_L^i + \pi_T^i, \quad \partial_i \pi_T^i = 0, \quad \epsilon_{ijk} \partial_j \pi_L^k = 0, \qquad (5.15)$$

where the transverse fields $\pi_T^i$ couple to vector metric perturbations and their evolution is dictated by the conservation of vorticity [10], which is a consequence of the symmetry $V_s$Diff. As we discuss in Appendix A, this can be easily obtained applying Noether's theorem, which implies and infinite set of covariantly conserved currents [9,10].

A peculiar property of the perfect fluids of the kind $U(b)$ is that the transverse modes $\pi_T^i$, which are interpreted as vortices, do not propagate in flat space since the symmetry $V_s$Diff forbids spatial derivatives in their action. Their time evolution is simply given by the equation $\ddot{\pi}_T^i = 0$. This has lead to the suggestion that these fluids may be strongly coupled at all scales [23]. However, it has been argued that the issues associated to the evolution of these modes (e.g. the appearance of divergences in the scattering of phonons in flat spacetime) can be resolved with a careful choice of the appropriate (physical) observables [32].

It can be easily checked, using a field redefinition, that the inclusion of the NLO operators respecting $V_s$Diff and

containing only the three spatial Stückelberg fields does not modify the evolution equation of $\pi_T^i$ in Minkowski (with respect to the one at LO) [10]. There are five independent operators at NLO: $\partial_\mu b \partial^\mu b$, $H^{\mu\nu} \nabla_\mu u^\alpha \nabla_\nu u_\alpha$, $(u^\mu \partial_\mu b)^2$, $\nabla_\mu u^\nu \nabla_\nu u^\mu$ and $\epsilon_{\alpha\beta\mu\nu} \nabla^\mu u^\alpha \nabla^\nu u^\beta$, see [10]. Once these operators are included, the EMT of the system is no longer that of a perfect fluid.

The relativistic formulation of self-gravitating media that we are exploring in this paper allows a straightforward connection to the thermodynamic theory of continuous media. This can be done by constructing a dictionary between two EFT pictures: the action (5.1) and the fundamental equations of thermodynamics, see e.g. [24]. A more detailed discussion of this dictionary is beyond the scope of this work and we provide it, for the case of perfect fluids, in [52]. For the purpose of illustration we use here the perfect fluid $U(b)$. Essentially, we need to match the operator $b$ with an intensive thermodynamic variable (or a combination of them). Choosing the particle number density $n$ and the entropy density $s$ as independent variables, we write $b = b(n, s)$. Then, the chemical potential $\mu$ and the temperature are obtained from the first principle of thermodynamics as $\mu = \partial\rho/\partial n$ and $T = \partial\rho/\partial s$. The pressure enters through the Euler relation

$$\mu s + nT = \rho + p, \qquad (5.16)$$

which must be valid for any $U(b)$. This can be achieved, for instance, with the identification $b = n$ (and then $T = 0$), which implies $\mu = -U_b$, and is in principle suitable for a degenerate system (i.e. in the limit $T \to 0$). A different option is $b = s$ (with $\mu = 0$) and $T = -U_b$, as for a gas of photons. This second choice is the one that was advocated e.g. in [10,23,24].

### 3. $U(X)$: Perfect fluid

At the opposite side of the spectrum from solids and perfect fluids lie the self-gravitating media that can be described at LO by introducing only the temporal Stückelberg field $\Phi^0$. We have not found a symmetry that forbids the spatial Stückelberg fields at LO but reintroduces them at NLO or at higher orders in derivatives. This is analogous to what happens for the solids and perfect fluids studied above, but exchanging the roles of the temporal and spatial fields. By looking at Table II, it seems that $V_s$Diff and $\Phi^a \to \Phi^a + f^a(\Phi^0)$ select $X$, but this combination of transformations is not a symmetry of an action based on $U(X)$, simply because $\Phi^a$ are not in $X$. Therefore, the absence of these fields must be assumed.

Clearly, with such an assumption, the function $U(X)$ is the only shift-symmetric possibility that exists at LO (with a single scalar field). Further restrictions involve necessarily discrete symmetries such as e.g. $X \to -X$.
Then for a generic $U(X)$, the EMT is given by

$$ 
\begin{aligned}
T_{\mu\nu} &= pg_{\mu\nu} + (\rho + p)\mathcal{V}_{\mu}\mathcal{V}_{\nu}, \\
\rho &= -U + 2XU_X, \quad p = U.
\end{aligned} 
\quad (5.17) 
$$

This type of perfect fluid is fundamentally different from the case $U(b)$ because its four-velocity, $\mathcal{V}^{\mu}$, is the derivative of a scalar and thus it has zero vorticity. This is due to the fact that for $U(X)$ the transverse phonons $\pi_T^i$ are entirely absent by construction. Actually, since both $U(X)$ and $U(b)$ describe perfect fluids, their phonon expansions (2.8) can be matched setting to zero the transverse modes of the latter and identifying through an equation the divergence of $\pi^i$ with the time derivative of $\pi^0$. This matching generically involves as well a linear combination of the scalar metric perturbation. The correspondence can also be easily established at the level of the background. For instance, a constant barotropic equation of state is obtained choosing $U \propto X^{(1+w)/(2w)}$, to be compared with (5.13). For an early study of irrotational perfect fluids where longitudinal fluctuations were already interpreted as phonons and a thermodynamic dictionary was provided, see [53].

In order to extend beyond LO the medium that can be constructed using only the temporal Stückelberg, we must linearly add to the action the operators $(\nabla_{\mu})\mathcal{V}^{\mu})^2$ and $\nabla_{\mu}\mathcal{V}^{\nu}\nabla^{\mu}\mathcal{V}_{\nu}$ at NLO. Clearly, the operator $\nabla^2\Phi^0 \sim \nabla_{\mu}\mathcal{V}^{\mu}$ can be omitted assuming appropriate boundary conditions, since it is a total derivative. Besides, the NLO operators $\epsilon_{\mu\nu\alpha\beta}\nabla^{\alpha}\mathcal{V}^{\mu}\nabla^{\beta}\mathcal{V}^{\nu}$ and $\mathcal{V}^{\mu}\mathcal{V}^{\nu}\nabla_{\mu}\mathcal{V}_{\alpha}\nabla_{\nu}\mathcal{V}^{\alpha}$ need not be included thanks to the fact that $\mathcal{V}^{\mu}$ is hypersurface orthogonal. The first of the two vanishes and the second one can be written as a combination of $(\nabla_{\mu}\mathcal{V}^{\mu})^2$ and $\nabla_{\mu}\mathcal{V}^{\nu}\nabla^{\mu}\mathcal{V}_{\nu}$.

If we assume that the operator $X$ is absent, and we work only at NLO, the resulting model is a special case of the “Einstein-aether” model [54], see e.g. [55].

Notice that the extension to all orders of $U(X)$ is not a generalization of the ghost condensate [56], but instead the most general modified gravity model based on a single scalar derivatively coupled. Actually, including all the orders in the derivative expansion corresponds to a large class of models contained in the EFT of inflation [57] (or DE [14]) for which the (soft) breaking of the temporal shift symmetry $\Phi^0 \to \Phi^0 + c^0$ (with constant $c^0$) due to a potential term for the inflation (or, generically, the extra scalar mode) is entirely neglected.

## B. Four-dimensional media

We will now move on to describe the media of Table II, which pertain to a different class than the cases we have studied so far, because all of them contain the four Stückelberg fields $\Phi^A$, which motivates the name we give them. As we will now see, the most general of these media (whose action is invariant under shifts and internal spatial rotations) is the covariant generalization of a nonrelativistic

model used in [46] to introduce supersolids. If we further constrain these systems by requiring invariance under internal spatial diffeomorphisms that preserve the volume, we get superfluids [24]. Symmetries that mix the temporal and spatial Stückelberg fields—see Table II—lead to a set of particular subcases of the general supersolid.

### 1. $U(X, Y, \tau_n, y_n)$: Supersolid

Let us start with the most general media that our symmetries allow. The less symmetric self-gravitating media that we can construct with our assumptions respect only the internal shifts (2.1) and the spatial rotations $SO(3)_s$ of (4.1). Therefore, their LO master function $U$ contains all the operators of the set (4.6), which we identified in Sec. IV, and thus their diffeomorphism invariant action is precisely (5.1). A nonrelativistic description of systems characterized precisely by these symmetries was given in [46], where they were interpreted as zero-temperature supersolids. According to [46], the three spatial $\Phi^a$ would correspond to the comoving coordinates of the medium (as we have been interpreting them), whereas $\Phi^0$ would be the (shift-symmetric) phase field related to a $U(1)$ symmetry associated to particle number conservation.

In spite of some claims, see e.g. [58], whether actual supersolids exist in nature is still unclear. For the purpose of this work we will simply use the term to refer to a continuous medium amenable to a coarse-grained description with four degrees of freedom and based on the shift and rotational symmetries (2.1) and (4.1). For a review on supersolids we point the reader to [59].

The same set of symmetries, Eqs. (2.1) and (4.1), was later considered in [30], where the CCWZ method [60] was applied to derive (for Minkowski spacetime) a relativistic version of the LO action given in [46]. However, the action written in [30] for these symmetries depends on a generic function of all our kinetic blocks $C^{AB}$ except $X$, which is an independent invariant operator under internal translations and $SO(3)_s$ rotations. Besides, the concrete operators $y_n$, defined in (4.4) were not given in [30].

Phenomenologically, the media of this kind possess many interesting properties. First of all, they share with the solids $U(\tau_n)$ the fact that the graviton acquires a nonvanishing mass that is related to the presence of anisotropic stress and to the speed of propagation of transverse modes. We refer to the group velocity $d\omega/dk$, which given e.g. a dispersion relation of the form $\omega^2 = k^2 + m^2$ clearly involves the mass $m$. In addition, given that $\Phi^0$ and $\Phi^a$ lead to two different four-vectors in the medium, the four-velocities (3.4) and (3.8), the energy-momentum tensor exhibits a nonvanishing “heat-flux”.

We point out that the object $C^{AB}$, defined in (2.2), is the four-dimensional generalization of the left Cauchy-Green tensor, see Sec. VA 1. This agrees with the picture in which the general media described by the EFT (5.1) can be
interpreted as four-dimensional media (or hypervolumes) propagating in spacetime. In Sec. VI we will make the connection between these media and massive (and modified) gravity, using the unitary gauge as tool to simplify the matching.

## 2. U(X,Y,b): Superfluid

Imposing that the action (5.1) has to be invariant under $V_sDiff$, the internal volume preserving diffeomorphisms (4.9), the only invariant LO operators that appear in the master function are $X$, $Y$ and $b$. The corresponding EMT is given by

$$T_{\mu\nu} = (U - bU_b)g_{\mu\nu} + (YU_Y - bU_b)u_{\mu}u_{\nu} + 2XU_XV_{\mu}V_{\nu}. \qquad (5.18)$$

In general this is not a perfect fluid, due to the four-velocities $V^\mu$ and $u^\mu$ not being parallel. The fact that $V^\mu$ has zero vorticity has motivated interpreting these media as relativistic superfluids [24,26], following earlier ideas for (nonrelativistic) supersolids [46] and superfluids [61]. The key idea of this picture is that the EMT (5.18) can be seen as a two-component fluid whose admixture allows us to describe the various phases of a superfluid. For instance, it is straightforward to show that the phonon expansion of the action leads to two kinds of longitudinal waves (propagating with different speeds of sound), which suggests interpreting them as the first and second sound in superfluids [26].

It is worth highlighting some special cases of the general superfluid.

(i) If in addition to the symmetry $V_sDiff$ we impose also invariance under $\Phi^0 \to \Phi^0 + f(\Phi^0)$, the master function has to be a function of the operators $O_a$, which we defined in (4.17), and $b$. The operators $O_{-1}$ and $b$, together with the NLO operators that can be constructed from $V^\mu$, were used in [41] to model Lorentz breaking in dark matter.

(ii) The master function $U(b, Y)$ describes a perfect fluid with four-velocity $u^\mu$ and energy density and pressure

$$\rho = YU_Y - U, \qquad p = U - bU_b. \qquad (5.19)$$

This case is interesting because it comes from imposing not only the symmetry $V_sDiff$ but also requiring invariance under $\Phi^0 \to \Phi^0 + f(\Phi^a)$, see Table II. This feature makes it qualitatively different from the other perfect fluids that we have encountered, $U(X)$ and $U(b)$, whose LO structure does not derive from symmetries acting on the four Stückelberg fields. One can obtain a constant barotropic equation of state for $U(Y, b)$ choosing

$$U \propto b^{1+w}U(b^{-w}Y), \qquad (5.20)$$

where $U$ is an arbitrary function. Another choice of master function which also leads to a constant barotropic equation of state is $U \propto (b^{1+w} + Y^{1+1/w})$.

(iii) An equation of state $w = -1$ can be obtained from (5.20) by choosing

$$U = U(bY), \qquad (5.21)$$

which interestingly corresponds to the enhanced symmetry

$$\Phi^A \to \Psi^A(\Phi^B), \qquad \det|\partial\Psi^A/\partial\Phi^B| = 1. \qquad (5.22)$$

By looking at Eq. (4.8), it is now clear why $bY$ is the operator that allows us to switch from Eulerian to Lagrangian coordinates.

## 3. Special supersolids

These media are the subclasses of the general action (5.1) listed in Table II.

(i) For instance, imposing the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^a)$, the leading invariant operators are $\tau_n$ and $Y$ and the EMT is

$$T_{\mu\nu} = U g_{\mu\nu} - 2U_{ab}\partial_{\mu}\Phi^a\partial_{\nu}\Phi^b + YU_Yu_{\mu}u_{\nu}. \qquad (5.23)$$

Given what we have seen so far, a natural interpretation of (5.23) seems to be that it describes a solid coupled to a perfect fluid. An analysis of the thermodynamics of perfect fluids [52] suggests that $Y$ can be interpreted as the temperature so that $U(\tau_n, Y)$ may be a possible description for a solid at finite temperature.

(ii) The symmetry $\Phi^a \to \Phi^a + f^a(\Phi^0)$ enforces the master function to be of the form $U = U(w_n, X)$. If in addition we impose the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^0)$, the operators $w_n$ are selected. We will comment on these cases in the next section, where we build the connection between massive gravity and self-gravitating media.

# VI. SO(3) SPATIALLY SYMMETRIC MASSIVE GRAVITY

The main motivation for massive and modified gravity in the context of cosmology is the possibility that the observed acceleration of the Universe could be due to a modification of GR that weakens gravity at large distances. In this respect, the idea of MG is appealing, since endowing the graviton with a mass could effectively make gravity a short range interaction (on cosmological scales).

The first attempt to provide a mass to the graviton dates back to 1939, when Fierz and Pauli [20] studied the most
general Lorentz invariant mass term for metric perturbations around Minkowski space. It was thought for a long time that any extension of the Fierz-Pauli Lagrangian beyond the quadratic level was pathological. Indeed, Boulware and Deser argued that at the nonlinear level a ghost mode would have necessarily be present [62]. The problem was solved in [63] by finding a ghost-free extension of the Fierz-Pauli Lagrangian at the nonlinear level [63,64]. This model, which exhibits Lorentz symmetry when expanded around Minkowski, solves the ghost issue at the (heavy) price of lacking spatially flat homogeneous FLRW solutions [65]. As shown in [66] the Fierz-Pauli model as an effective theory (and also [63]) has a (very low) cutoff $\Lambda_3 = (m^2 M_{pl})^{1/3}$, where m sets the scale of the graviton mass scale. A recent proposal to improve this cutoff has been given in [67,68].

Fierz-Pauli MG fails to reproduce the correct light bending around heavy and dense objects, in sharp contrast with standard GR. This also happens in the limit of vanishing graviton mass, an issue that is known as the vDVZ discontinuity [69]. A possible way out was proposed by Vainshtein in [70], by arguing that nonlinear effects restore the correct GR behavior at short distances from a gravitational source. Whereas this mechanism has been studied for various models [71], it is important to stress that it relies on strong nonlinearities, even at the macroscopic scales of the Solar System where the value of the gravitational potential is small.

A possible way out of this last difficulty is to abandon Lorentz invariance in the gravitational sector<fnref n="9" /> in favor of a simpler symmetry group, in practice, by requiring only spatial rotational invariance [7,73]. All Lorentz-breaking MG models with five degrees of freedom (d.o.f.) having (at least) spatial rotational invariance are free of ghost instabilities and can be classified [11,12]. In addition, they have no vDVZ discontinuity and FLRW solutions do exist for these models [74]. Their range of validity is given by an ultraviolet cutoff $\Lambda$, which is dictated by the symmetries that define the EFT. Typically, for Lorentz-breaking MG the cutoff is $\Lambda_2 = (mM_{pl})^{1/2} \gg \Lambda_3$ [7,8,73,74]. In order to have any impact on the current acceleration of the Universe, the scale m should thus be of the order of today's Hubble scale, giving $\Lambda_2^{-1} \sim 0.1$ mm.

The EFT of four derivatively coupled scalar fields with internal spatial $SO(3)_s$ invariance allows us to reinterpret Lorentz-breaking MG models as self-gravitating media. From the perspective of this EFT, the breaking of Lorentz symmetry is simply a consequence of a background choice. MG is then a class of models within a broad framework. In other words, the EFT of self-gravitating media leads to a

systematic construction of models for the acceleration of the Universe, including MG.

The basic idea that helps to build this connection is the Stückelberg mechanism or trick, which we already mentioned in Sec. II. The application of the Stückelberg mechanism for Fierz-Pauli MG was already given in [75] and was later generalized to other models of MG in [66] and [7]. If $h_{\mu\nu}$ is a metric perturbation around a (reference) background metric, $\tilde{g}_{\mu\nu}$, such that the full metric is $g_{\mu\nu} = \tilde{g}_{\mu\nu} + h_{\mu\nu}$, a generic (Lorentz invariant) MG potential may be expressed as function of the perturbation $h_{\mu\nu}$ with indices raised using $g^{\mu\nu}$, i.e. $\sqrt{-\tilde{g}}V(h_{\mu\nu}, g^{\mu\nu})$. Given such a potential, added to the Einstein-Hilbert Lagrangian, the Stückelberg trick can be implemented replacing the background metric $\tilde{g}_{\mu\nu}$ with a tensor field made out of four scalars, $G_{\mu\nu} \equiv \partial_{\mu}\Phi^A \partial_{\nu}\Phi^B \tilde{g}_{AB}[\Phi^C(x^\alpha)]$, so that the metric perturbation $h_{\mu\nu}$ is replaced with $g_{\mu\nu} - G_{\mu\nu}$, see e.g. [76]. Using these replacements to express the action in terms of the spacetime metric $g_{\mu\nu}$ and the tensor $G_{\mu\nu}$ we obtain a covariant embedding of any Lorentz invariant MG model originally defined with a reference metric $\tilde{g}_{\mu\nu}$. This procedure can be adapted to deal also with Lorentz breaking models, as described in [7].

In this work, we started our analysis directly with an explicitly diffeomorphism invariant action (2.3) in a “top-down” approach [8]. This automatically implements the Stückelberg mechanism for the actions obtained imposing the background $\Phi^A = \delta^A_\mu x^\mu$. Indeed, assuming that the condition $\det(\partial_\mu \Phi^A) \neq 0$ holds, it is possible to find a local set of coordinates (unitary gauge) such that

$$ \partial_{\mu} \Phi^{A} = \delta_{\mu}^{A}. \qquad (6.1) $$

Notice that, defined in this way, the unitary gauge is actually a collection of gauges, which are equivalent modulo constant shifts, i.e. $\Phi^A = x^A + c^A$ with $\partial c^A / \partial x^\mu = 0$, where for simplicity we do not distinguish between spacetime and internal manifold indices. These gauges are all equivalent since the action is, by assumption, shift invariant in the Stückelberg fields.

The unitary gauge reflects the fact that since we have four scalars $\Phi^A$ and four spacetime coordinates, we can use the diffeomorphism invariance of the action (5.1) to absorb the dynamics of the scalar fields into the degrees of freedom of the metric. By using the ADM splitting of spacetime to describe the metric $g_{\mu\nu}$ and its perturbations [77], it is then straightforward to make a connection with MG in the so-called broken phase. Therefore, this phase (of broken diffeomorphisms) can be simply understood as a way of expressing the dynamics of the comoving coordinates of the self-gravitating medium through their effect on the metric perturbations. The correspondence between self-gravitating media and massive and modified gravity models written explicitly in terms of metric perturbations

<footnote n="9">⁹Bounds on Lorentz violation in the gravitational sector are rather mild assuming the equivalence principle and come from post-Newtonian preferred frame effects, see e.g. [72], and gravitational wave emission.</footnote>
<table><caption>TABLE IV. Relationship between material Lagrangians and massive gravity models.</caption><thead><tr><th>LO self-gravitating media</th><th>Map</th><th>Massive gravity</th></tr></thead><tbody><tr><td><i>L</i>(C<sup>AB</sup>, g<sub>μν</sub>)</td><td>Unitary gauge →</td><td><i>L</i>(h<sub>μν</sub>, g<sup>μν</sup>)</td></tr><tr><td><i>O</i><sub>LO</sub>: X, Y, τ<sub>n</sub>, y<sub>n</sub></td><td>← Stückelberg “trick”</td><td>SO(3) invariants of ADMs N, N<sup>i</sup>, γ<sub>ij</sub></td></tr></tbody></table>

is depicted in Table IV, where the arrows indicate how to move from one picture to the other.

In the ADM formalism the metric can be written in terms of the lapse N, the shifts N<sup>i</sup> and the spatial metric γ<sub>ij</sub>, namely

$$ 
\begin{aligned}
g_{\mu\nu} &= \begin{pmatrix} -N^2 + N^i N^j \gamma_{ij} & \gamma_{ij} N^j \\ \gamma_{ij} N^j & \gamma_{ij} \end{pmatrix}, \\
g^{\mu\nu} &= \begin{pmatrix} -1/N^2 & N^i/N^2 \\ N^j/N^2 & \gamma^{ij} - N^i N^j/N^2 \end{pmatrix},
\end{aligned} 
\quad (6.2) 
$$

where γ<sup>im</sup>γ<sub>mj</sub> = δ<sup>i</sup><sub>j</sub>. In the unitary gauge (6.1), the matrix C<sup>AB</sup> defined in (2.2) coincides with inverse of the spacetime metric,

$$ C^{AB} = g^{\mu\nu} \delta_{\mu}^{A} \delta_{\nu}^{B}. \quad (6.3) $$

Defining

$$ \xi^i \equiv N^i/N \quad (6.4) $$

for convenience, we decompose C<sup>AB</sup> in the unitary gauge as follows:

$$ X = C^{00} = -1/N^2, \quad C^i = C^{i0} = \xi^i/N, \quad B^{ij} = C^{ij} = \gamma^{ij} - \xi^i \xi^j. \quad (6.5) $$

Then, we can write the building blocks of LO self gravitating media—see (4.6) and Table I—in terms of the ADM variables,

$$ 
\begin{aligned}
&X = -1/N^2, \qquad Y = \frac{1}{N\sqrt{1 - \xi^2}}, \\
&\tau_n = \mathrm{Tr}[(\gamma - \xi \otimes \xi)^n], \qquad y_n = \frac{\xi \cdot (\gamma - \xi \otimes \xi)^n \cdot \xi}{N^2},
\end{aligned} 
\quad (6.6) 
$$

where γ is the 3 × 3 matrix of components γ<sup>ij</sup>. In these expressions the dot (·) represents the standard matrix product and we use the following notation: ξ<sup>2</sup> ≡ ξ<sup>i</sup>γ<sub>ij</sub>ξ<sup>j</sup> and (ξ ⊗ ξ)<sup>ij</sup> ≡ ξ<sup>i</sup>ξ<sup>j</sup>. Therefore, the master function U of the action (5.1) becomes an algebraic function of the ADM variables in the unitary gauge,

$$ S = M_{pl}^2 \int d^4 x \sqrt{-g} R + \int d^4 x \sqrt{-g} U(N, \xi^i, \gamma^{ij}). \quad (6.7) $$

Special cases—see Tables II and III—can be obtained from the remaining SO(3)<sub>s</sub> invariant scalar operators of Table I, which are combinations of X, Y, τ<sub>n</sub> and y<sub>n</sub>,

$$ 
\begin{aligned}
&b = \sqrt{(1 - \xi^2) \det \gamma}, \\
&w_n = \mathrm{Tr}(\gamma^n), \quad \mathcal{O}_{\alpha\beta n} = (-1)^\alpha (1 - \xi^2)^{\alpha+\beta} (\xi \cdot (\gamma - \xi \otimes \xi)^n \cdot \xi)^\beta.
\end{aligned} 
\quad (6.8) 
$$

Notice also that the four-velocities V<sup>μ</sup> and u<sup>μ</sup> become in the unitary gauge

$$ u^{\mu} = \frac{\delta_{0}^{\mu}}{N\sqrt{1 - \xi^{2}}}, \qquad \mathcal{V}^{\mu} = (-1/N, \xi^{i}). \quad (6.9) $$

If the shifts are zero, e.g. as in Minkowski and FLRW spacetimes, u<sup>μ</sup> and V<sup>μ</sup> coincide, up to a sign; see also Sec. III.

The action (6.7) was the starting point in [11–13] where a large class of nonderivative<sup>10</sup> massive and modified gravity models—those with spatial SO(3) invariance—was studied by means of a Hamiltonian analysis. Models with six, five, three and two (d.o.f.) were thus found. The condition of SO(3) spatial invariance turns to be important to avoid ghost instabilities. If this restriction is relaxed, six d.o.f. generically propagate and one of them (a scalar mode with respect to spatial rotations) is a ghost around flat space, which corresponds to the infamous Boulware-Deser ghost [62].

For example, a SO(3) invariant master function U of the form

$$ U = U(\tau_n) + \sqrt{-X} \mathcal{E}(w_n) \quad (6.10) $$

propagates precisely five d.o.f. [12] (provided that U is not a constant). As a matter of fact, this is not the most general U that propagates five d.o.f., see [12]. We leave open a detailed exploration of the structure of the most general case, which is contained in our formalism, for a possible future investigation. The functions U and E of (6.10) are

<footnote n="10"><sup>10</sup>By “nonderivative” we refer to the actions for which the master function U contains only nonderivative combinations of the ADM variables, precisely as in (6.7).</footnote>
generic and thus there exist infinitely many U with five d.o.f. Therefore, according to the discussion in the previous section ghost-free massive gravity can be identified with a specific kind of supersolid. Moreover, setting $\mathcal{E} = 0$ there are still five d.o.f. and the medium is a solid. Thus, the construction of ghost-free massive gravity with five d.o.f. needs at least the three spatial Stückelberg fields $\Phi^a$.

Cases where only three or two d.o.f. occur also exist [13]. In particular, exactly three d.o.f. are present if U has the structure

$$U = U(w_n, X). \qquad (6.11)$$

We recall that the above form for U is protected by the symmetry $\Phi^a \to \Phi^a + f^a(\Phi^0)$, see Table II. It should be stressed that although three d.o.f. are present at the non-perturbative level for (6.11), only two d.o.f. are found expanding around flat space [13]. A proposal for a UV completion, where the LO order part of the action is built out of the scalar operators $w_1$ and $w_2$, was put forward in [78] adding the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^0)$, see Table II. In that model the fields $\Phi^a$ are further coupled with a triplet of purely spatial vector fields, whereas the dynamics $\Phi^0$ is dictated by NLO operators.

A particular case of (6.11) is given by the Lagrangian

$$U = \sqrt{-X}\mathcal{E}(w_n) + \lambda, \qquad (6.12)$$

where only two nonperturbative d.o.f. are found. In this expression $\lambda$ is a (dimensionless) cosmological constant. This last case is rather interesting since it gives an example of a model of gravity with two d.o.f. that is different from GR.

## VII. COSMOLOGICAL PERTURBATIONS IN FLRW

In this section we discuss a basic aspect of cosmological perturbation theory for self-gravitating media. In particular, we are interested in determining the number of degrees of freedom that propagate in a FLRW background for the different kinds of self-gravitating media that we have considered. A more detailed account of the cosmology of these media will be given elsewhere [79]. Here we present the Lagrangians for perturbations that are required for such a study and discuss their main features.

Let us consider a spatially flat<fnref n="11" /> FLRW metric written in the following way:

$$ds^2 = -N(t)^2 dt^2 + a(t)^2 \delta_{ij} dx^i dx^j, \qquad (7.1)$$

where we denote $x^0 = t$. The derivatives with respect to $t$ will be indicated by primes. For instance, by definition

$\mathcal{H} = a'/a$ is the Hubble function in these coordinates. With this choice of metric, the field configuration

$$\Phi^0 = \phi(t), \qquad \Phi^i = x^i \qquad (7.2)$$

is compatible with the equations of motion for the Stückelberg fields derived from the action (5.1), giving a spatially homogeneous and isotropic EMT. In what follows we will work in the unitary gauge, setting $\phi(t) = t$. It is important to stress that in general it is not possible to choose $N(t) = a(t)$ or $N(t) = 1$ at the same time that $\phi(t) = t$. Therefore, $N(t)$ has to be determined using the (background) equations of motion.

We will now study perturbations around the metric (7.1). Choosing the unitary gauge, the dynamics of perturbations for the whole system (metric and Stückelberg fields) is encoded in the metric.<fnref n="12" /> Denoting the metric perturbations by $h_{\mu\nu}$, the master function U of the action (5.1) can be expanded up to second order in $h_{\mu\nu}$ as follows:

$$\sqrt{-g} U = t^{\mu\nu} h_{\mu\nu} + \frac{M_{pl}^2}{4} (m_0^2 h_{00}^2 + 2m_1^2 h_{0i} h_{0i} - 2m_4^2 h_{00} h_{ii} + m_3^2 h_{ii}^2 - m_2^2 h_{ij} h_{ij}), \qquad (7.3)$$

where $t^{\mu\nu}$ and the masses $m_i^2$ depend on U and its derivatives.

Although the operator b is redundant once the $\tau_n$ are included, it is convenient to single out its effect, since it plays a special role for some symmetry choices. Hence, we write the LO action (5.1) as

$$S = M_{pl}^2 \int d^4 x \sqrt{-g} R + \int d^4 x \sqrt{-g} U(X, Y, \tau_n, y_n, b). \qquad (7.4)$$

Then the masses of (7.3) are given by

$$m_0^2 = -\frac{a^7}{2N^3 M_{pl}^2} \left[ \frac{4U_{XY}}{N^3} - \frac{4U_{XX}}{N^4} - \frac{U_{YY} - 4U_X}{N^2} - \frac{U_Y}{N} + U \right], \qquad (7.5)$$

$$m_1^2 = \frac{a}{NM_{pl}^2} \left[ a^4 U + \frac{2a^4}{N^2} U_X - aU_B + \sum_{n=0}^{3} 2na^{2-2n} U_{Y_n} - 2 \sum_{n=1}^{3} a^{4-2n} U_{\tau_n} \right], \qquad (7.6)$$

$$m_2^2 = \frac{N}{M_{pl}^2} \left[ a^3 U - U_B - 4 \sum_{n=1}^{3} (n+1) a^{3-2n} U_{\tau_n} \right], \qquad (7.7)$$

<footnote n="11">¹¹Spatial curvature can be introduced without altering the conclusions. For simplicity we set it to zero.</footnote>

<footnote n="12">¹²In the unitary gauge, the “phonons” $\pi^A$ of (2.8) are set to zero. We could as well choose to trade some degrees of freedom of the metric (7.15) by the $\pi^A$.</footnote>
$$m_3^2 = \frac{N}{M_{pl}^2} \left[ \frac{1}{2}(a^3 U + a^{-3} U_{bb} - U_b) + 8a^{-3}U_{\tau_1\tau_2} + 12a^{-5}U_{\tau_1\tau_3} + 24a^{-7}U_{\tau_2\tau_3} \right. \\ \left. + 2\sum_{n=1}^3 (na^{-2n}U_{b\tau_n} - na^{3-2n}U_{\tau_n} + n^2a^{3-4n}U_{\tau_n\tau_n}) \right], \qquad (7.8)$$

$$m_4^2 = \frac{a^2}{2N^3 M_{pl}^2} \left[ a^3 N^2 U - a^3 N U_Y + 2a^3 U_X - N^2 U_b - 2U_{bX} + N U_{bY} \right. \\ \left. + 2\sum_{n=1}^3 na^{3-2n}(NU_{Y\tau_n} - N^2 U_{\tau_n} - 2U_{X\tau_n}) \right]. \qquad (7.9)$$

These masses can be simplified (on a case by case basis) taking into account the background equations of motion. Concretely, the condition $t^{\mu\nu} = 0$ is required for consistency. This gives two equations,

$$6M_{pl}^2\mathcal{H}^2 + N^2U - NU_Y + 2U_X = 0,$$

$$N^3 U + 2M_{pl}^2 N(2\mathcal{H}' + 3\mathcal{H}^2) - 4M_{pl}^2 \mathcal{H}N' - 2N^3 \sum_{n=1}^{3} na^{-2n}U_{\tau_n} - a^{-3}N^3 U_b = 0. \qquad (7.10)$$

These equations are equivalent to the conservation of the medium's EMT on the background or, equivalently, to the equations of motion of the Stückelberg fields (in the background), see Sec. II. Eliminating $\mathcal{H}'$ from the second of them we obtain the condition

$$F_1(U, N, a)N' + \mathcal{H}F_2(U, N, a) = 0, \qquad (7.11)$$

where

$$F_1 = -\frac{a^3}{2N^4\mathcal{H}}(N^2(U_{YY} - 2U_X) - 4NU_{XY} + 4U_{XX}), \qquad (7.12)$$

$$F_2 = \frac{3}{2} \left( a^3 U_Y - U_{bY} - 2 \sum_{n=1}^{3} na^{3-2n} U_{\tau_n} \right) \\ + \frac{3}{N} \left( U_{bX} - aU_X + 2 \sum_{n=1}^{3} na^{3-2n} U_{X\tau_n} \right). \qquad (7.13)$$

Therefore, we have the following possibilities:

(i) $F_1 \neq 0$ and therefore it is possible to solve for $N'$. As a result, $N$ is dynamical.

(ii) $F_1$ is identically zero. Therefore, (7.11) is an algebraic equation for $N$, i.e. $F_2(U, N, a) = 0$, except if $\mathcal{H} = 0$.

(iii) Both $F_1$ and $F_2$ are identically zero. In this case $N$ is not fixed by the background (so there exists a residual gauge freedom). We can fix $N$, for instance, to be equal to $a$.

Let us now consider the scalar, vector and tensor modes of the metric with respect to spatial rotations. A general perturbation of (7.1) can be written as

$$ds^2 = -N(t)^2 dt^2 + a(t)^2 \delta_{ij} dx^i dx^j + a(t)^2 h_{\mu\nu} (x^\alpha) dx^\mu dx^\nu, \qquad (7.14)$$

where

$$h_{00} = \psi, \qquad h_{ij} = \chi_{ij} + \partial_i s_j + \partial_j s_i + \delta_{ij} \tau + \partial_i \partial_j \sigma, \\ h_{0i} = u_i + \partial_i v, \qquad (7.15)$$

with

$$\partial_i u_i = \partial_j \chi_{ij} = \partial_i s_i = 0, \qquad \delta_{ij} \chi_{ij} = 0. \qquad (7.16)$$

In the scalar sector, one can integrate out the fields $v$ and $\psi$ of (7.15), arriving at the following canonical form for the effective quadratic Lagrangian

$$L_{(2)}^{(s)} = \frac{1}{2} \varphi^{\prime t} K \varphi^{\prime} + \varphi^{\prime t} D \varphi - \frac{1}{2} \varphi^{\prime} A_{\text{mass}} \varphi, \qquad \varphi^{\prime} = (\tau, -\Delta\sigma), \qquad (7.17)$$

where the $2 \times 2$ matrices $K$, $D$ and $A_{\text{mass}}$ satisfy $D^t = -D$, $K^t = K$ and $A_{\text{mass}}^t = A_{\text{mass}}$, respectively. The propagation of modes is related to the determinant of the matrix $K$: if $\det(K) = 0$ there can be at most one propagating scalar d.o.f. We have that

$$\det(K) \propto (6a^5\mathcal{H}^2 + m_1^2 N^3)(3a^7\mathcal{H}^2 - m_0^2 N^5). \qquad (7.18)$$
Defining the effective masses

$$M_0^2 = m_0^2 N^5 - 3a^7 \mathcal{H}^2, \quad M_1^2 = 6a^5 \mathcal{H}^2 + m_1^2 N^3, \quad (7.19)$$

we see that if $M_0^2 = 0$ or $M_1^2 = 0$, there is a single propagating scalar mode, whereas if $M_1^2 = M_0^2 = 0$ no scalar mode propagates. In a general case, two different scalar modes can propagate. The study of scalar perturbations is important for cosmological signatures that could help to distinguish a modification of gravity such as the ones discussed here from a pure cosmological constant, see [80] for the prospects with the Euclid satellite. First, given a certain background evolution, the growth of dark matter perturbations is generically altered if there are extra scalar modes. This can be detected by measuring precisely the so-called growth index and its scale and redshift dependence, see e.g. [81–85]. Moreover, a detection of the speed of propagation of DE perturbations would also indicate a deviation from the $\Lambda$CDM model, see e.g. [86–91]. Besides, a late-time measurement of a non-negligible excess in the gravitational slip (i.e. the difference between the two gauge-invariant Bardeen potentials [92]) would signal the presence of an anisotropic stress component, see e.g. [83,93,94]. In our case, a gravitational slip appears in solids and supersolids, for which (at linear order in perturbations)

$$\delta T_{ij} \supset M_2^2 \partial_i \partial_j \pi_L, \quad (7.20)$$

where $\pi_L$ is defined via $\pi_L^i = \partial_i \pi_L$ and $\pi^i = \pi_L^i + \pi_T^i$, being $\epsilon_{ijk} \partial_j \pi_L^k = 0$, see (2.8), and we define

$$M_2^2 M_{pl}^2 = \sum_{n=1}^{3} n^2 a^{-2(n-1)} U_{\tau_n}. \quad (7.21)$$

Notice also that in a general case, the two scalar d.o.f. that appear from the action (7.4) allow for the propagation of an adiabatic mode and an entropy mode, see also [52].

In the tensor sector there are always two propagating d.o.f., with quadratic Lagrangian

$$L_{(2)}^{(t)} = \frac{M_{pl}^2}{2} \left[ \frac{a^3}{N} \chi'_{ij} \chi'_{ij} + Na(2M_2^2 - k^2) \chi_{ij} \chi_{ij} \right]. \quad (7.22)$$

A nonzero $M_2^2$ changes the propagation speed of gravitational waves with respect to that of light. This could induce, for instance, observable effects on the propagation and lensing of cosmic microwave background B-modes if the continuous medium is relevant at sufficiently early times, see [95,96]. It is remarkable that the same quantity, $M_2^2$, that is responsible for the gravitational slip also determines the propagation of gravitational waves. The relation between these two effects has also been noted for VDiff systems of three scalars at NLO in derivatives [16] and for other examples of modified gravity models [97].

Finally, in the vector sector the quadratic Lagrangian reads

$$L_{(2)}^{(v)} = M_{pl}^2 \left[ \frac{a^3 k^2}{2N} (u_i - s_i')(u_i - s_i') + \frac{M_1^2}{2N^3} u_i u_i + a N k^2 M_2^2 s_i s_i \right]. \quad (7.23)$$

The fields $u_i$ have a purely algebraic equation of motion and thus can be integrated out, giving the Lagrangian

$$\tilde{L}_{(2)}^{(v)} = k^2 M_{pl}^2 \left[ \frac{a^3 M_1^2}{2a^3 N^3 k^2 + 2N M_1^2} s_i' s_i' + Na M_2^2 s_i s_i \right]. \quad (7.24)$$

The vectors $s_i$ propagate only if $M_1^2 \neq 0$, see Table V. The dispersion relation is not trivial only when $M_2^2 \neq 0$. Thus, $M_2^2$, besides controlling the dispersion relation of tensors, also determines the dynamics of the vectors. The fact that $M_2^2 \neq 0$ in the unitary gauge is equivalent to the presence of an anisotropic stress in the EMT of the Stückelberg fields. As we mentioned above, a complete study of the cosmology of these media will be given in a separate publication [79]. We summarize in Table V the d.o.f. counting for various examples of sets of operators that may appear inside the LO master function.

<table><caption>TABLE V. Counting of d.o.f. for various LO self-gravitating media in FLRW.</caption><thead><tr><th>Media</th><th>Operators</th><th>Det(K)</th><th>M<sub>0</sub><sup>2</sup></th><th>M<sub>1</sub><sup>2</sup></th><th>M<sub>2</sub><sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td rowspan="4">Perfect fluids</td><td>b</td><td>0</td><td>0</td><td>&ne; 0</td><td>0</td><td>3</td></tr><tr><td>X</td><td>0</td><td>&ne; 0</td><td>0</td><td>0</td><td>3</td></tr><tr><td>Y</td><td>&ne; 0</td><td>&ne; 0</td><td>&ne; 0</td><td>0</td><td>3</td></tr><tr><td>b, Y</td><td>&ne; 0</td><td>&ne; 0</td><td>&ne; 0</td><td>0</td><td>3</td></tr><tr><td rowspan="2">Superfluids</td><td>b, X</td><td>&ne; 0</td><td>&ne; 0</td><td>&ne; 0</td><td>0</td><td>4</td></tr><tr><td>b, X, Y</td><td>&ne; 0</td><td>&ne; 0</td><td>&ne; 0</td><td>0</td><td>4</td></tr><tr><td>Solid</td><td>&tau;<sub>n</sub></td><td>0</td><td>0</td><td>&ne; 0</td><td>&ne; 0</td><td>5</td></tr><tr><td rowspan="2">Special supersolid</td><td>&tau;<sub>n</sub>, Y</td><td>&ne; 0</td><td>&ne; 0</td><td>&ne; 0</td><td>&ne; 0</td><td>6</td></tr><tr><td>w<sub>n</sub>, X</td><td>0</td><td>&ne; 0</td><td>0</td><td>&ne; 0</td><td>2</td></tr></tbody></table>
VIII. DARK MATTER AND DARK ENERGY SIGNATURES

To illustrate the novel phenomenology that arises for the description of dark energy and dark matter from the effective theory of self-gravitating media, we will focus now on certain aspects of cosmological perturbations that are specific of our framework. In particular, we will consider scalar perturbations for nonbarotropic fluids and also the propagation of gravitational waves.

A. Scalar perturbations in nonbarotropic fluids

As we mentioned in Sec. V B 2, if the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^a)$ is imposed on the general superfluid master function $U(b, X, Y)$, the operator $X$ drops and the resulting $U(b, Y)$ system is a perfect fluid. This type of perfect fluid is generically nonbarotropic because the pressure and the energy density cannot be written one as a function of the other, see (5.19). This kind of action can be used to describe a generic perfect fluid, since $b$ can be identified with the entropy density and $Y$ with the chemical potential.<fnref n="13" /> Then, the entropy per particle $\sigma = b/U_Y$ is nonlinearly conserved in time [52]. The pressure perturbation can be written as

$$ \delta p = c_s^2 \delta \rho + \frac{\phi'}{a^4} (c_s^2 - c_b^2) \delta \sigma, \quad c_s^2 = \frac{p'}{\rho'}, \quad c_b^2 = \left. \frac{\partial p}{\partial \rho} \right|_b, \qquad (8.1) $$

We recall that the quantity $\Gamma \equiv \delta p - c_s^2 \delta \rho$ (that is sometimes called “intrinsic entropy perturbation”) is gauge invariant and $\psi$ corresponds to one of the (gauge-invariant) Bardeen potentials [92]. Equation (8.4) is a linear combination of the 0-0 and the $i-i$ Einstein equations. In the case at hand, it contains a generic time-dependent speed of sound $c_s^2$, plus a time-dependent source term that is due to the nonbarotropic character of $U(b, Y)$.

Indeed, whereas $\delta\sigma \equiv \delta\sigma(k)$ is constant in time, with its $k$-dependence fixed by initial conditions, the factor that multiplies it in Eq. (8.4) has a nontrivial time dependence determined by $c_s^2$ and $c_b^2$, which in turn depend on the function $U(b, Y)$. Dynamical (and thermodynamic) stability of the perturbations requires $c_s^2 > 0$; see [52].

where $\phi = \phi(t)$ describes the background evolution of the field $\Phi^0$, see (7.2). In this expression and in all that follows, primes denote derivatives with respect conformal time $t$. The quantities $c_s^2$ and $c_b^2$ are, respectively, the standard adiabatic sound speed (i.e. the variation of the pressure with respect to the density at fixed $\sigma$) and the variation of the pressure with respect to the density at constant $b$. The latter of the two determines the evolution of $\phi$, which obeys (from the conservation of the EMT)

$$ \phi'' + \mathcal{H}(3c_b^2 - 1)\phi' = 0. \qquad (8.2) $$

Let us now write the scalar part of the perturbed metric in the conformal Newtonian gauge,

$$ ds^2 = -a^2(t)(1-\psi)dt^2 + a^2(t)(1+\tau)\delta_{ij}dx^i dx^j. \qquad (8.3) $$

We will work in Fourier space, so that the perturbations $\psi$ and $\tau$ will be functions of conformal time and the modulus, $k$, of the comoving momentum $k^i$. $U(b, Y)$ being a perfect fluid, there is no anisotropic stress and, in consequence, $\tau = \psi$. The evolution of the metric perturbation is then governed by the equation

$$ \psi'' + 3\mathcal{H}(c_s^2 + 1)\psi' + [c_s^2 k^2 + (3c_s^2 + 1)\mathcal{H}^2 + 2\mathcal{H}']\psi - \frac{\phi'(c_s^2 - c_b^2)}{2a^2 M_{pl}^2} \delta\sigma = 0. \qquad (8.4) $$

In the special case of nonzero and constant $c_b^2$, from (8.2) one has that

$$ \phi' = \phi_0' a^{1-3c_b^2}, \qquad (8.5) $$

where $\phi_0'$ is a constant and we take $a = a_0 = 1$ today. Then, the solution of (8.4) during a matter domination period, where $\mathcal{H} \propto 1/\sqrt{a}$, is given by

$$ \psi = \tilde{\psi} - \frac{\phi_0' \delta \sigma}{3 M_{pl}^2 H_0^2 \Omega_m^0 (6c_b^2 - 5)c_b^2} a^{-3c_b^2}, \qquad (8.6) $$

where $\tilde{\psi} = \tilde{\psi}(k)$ is the constant solution of the homogeneous equation and $\Omega_m^0$ and $H_0 = \mathcal{H}_0$ are, respectively, the relative matter density and the Hubble parameter in cosmic time, both evaluated today. We have neglected the decreasing solution of (8.4) and the effect of $c_s^2$, which is assumed to be positive but sufficiently small to be negligible for dark matter.

The energy density perturbation can be written directly from the 0-0 Einstein equation,

<footnote n="13">¹³This thermodynamic interpretation is not unique: $b$ can alternatively be associated to particle number density, and $Y$ to temperature, see [52]. However, we stress that the results of this section are independent of this choice and can be derived purely algebraically, without resorting to either of the two possibilities.</footnote>
$$ \delta\rho = 2a^{-2}M_{pl}^2[\psi(k^2 + 3H^2) + 3H\psi'] . \quad (8.7) $$

We get from it the time evolution of the density contrast
$$ \delta = \delta\rho/\rho, $$

$$ \delta = \delta_{ad} - \frac{\hat{\sigma}}{3(6c_b^2 - 5)c_b^2} \left( 1 - 3c_b^2 + \frac{ak^2}{3H_0^2\Omega_m^0} \right) a^{-3c_b^2}, \quad (8.8) $$

where

$$ \delta_{ad} = \tilde{\psi} \left( 1 + \frac{ak^2}{3H_0^2\Omega_m^0} \right) \quad (8.9) $$

is the nonentropic (or adiabatic) part of the fluctuation and we define the constant

$$ \hat{\sigma} = \frac{\phi_0' \delta\sigma}{M_{pl}^2 H_0^2 \Omega_m^0}. \quad (8.10) $$

With respect to the standard cold dark matter (for which $\delta\sigma = 0$), the entropic term modifies the time evolution of the density contrast. It is interesting to note that as soon as $c_b^2 < 0$, there is an entropic mode that grows with respect to the adiabatic one, both for super- and subhorizon scales.<fnref n="14" /> As an explicit example we can take the function $U(b, Y)$ corresponding to a perfect nonrelativistic gas [52],

$$ U(b, Y) = bY \left\{ 1 + \mathrm{Log} \left[ \frac{g_s}{b} \left( \frac{mY}{2\pi} \right)^{3/2} \right] \right\} - bm, \quad (8.11) $$

where $g_s$ is the number of spin states and $m$ is the atomic mass of the gas. One gets $c_b^2 = 2/3$ so that $\phi' \propto a^{-1}$, see (8.5), and

$$ \delta = \delta_{ad} + \frac{\hat{\sigma}}{2} \left( -1 + \frac{ak^2}{3H_0^2\Omega_m^0} \right) a^{-2}. \quad (8.12) $$

Therefore, in this example the entropic corrections decrease in time with respect the adiabatic perturbation (since $c_b^2 > 0$).

A more dramatic example of the influence of $\delta\sigma$ is obtained, for instance, by taking

$$ U(b, Y) = -\lambda b + V(bY). \quad (8.13) $$

In this case, $c_b^2 = -1$ and $\phi' \propto a^4$, see (8.5), so that $bY = \phi_0'$ is constant in time, and

$$ c_s^2 = 0, \quad w = p/\rho = \frac{V(\phi_0') - \phi_0' \frac{dV}{d(bY)}|_{\phi_0'}}{\lambda/a^3 - V(\phi_0') + \phi_0' \frac{dV}{d(bY)}|_{\phi_0'}}. \quad (8.14) $$

In the regime in which $\lambda/a^3$ dominates the denominator, $w \to 0$ (as cold dark matter) in a first approximation. In the opposite limit, $w \to -1$ and the fluid describes a phase of dark energy domination. During matter domination, when $w \sim 0$, the energy density contrast is given by

$$ \delta = \delta_{ad} - \frac{\hat{\sigma}}{33} \left( 4 + \frac{ak^2}{3H_0^2\Omega_m^0} \right) a^3. \quad (8.15) $$

In this case, the entropic perturbation grows faster than the adiabatic one. Although in this example, in order to keep under control the validity of linear perturbation theory, very strong constraints have to be imposed on the initial conditions (specifically on $\hat{\sigma}$), it illustrates the relevance of possible interactions between dark matter and dark energy in the context of the effective theory of self-gravitating media.

In general, by splitting $U(b, Y)$ as $U_{DM}(b, Y) + V(bY)$ one can systematically study the coupling of dark matter and dark energy, the corresponding EMT tensors not being separately conserved.

To conclude this subsection it is also worth hinting how the evolution of perturbations is affected for a superfluid and a solid. In the first of these two cases, described by a function $U(b, X, Y)$, the function $\delta\sigma$ is not constant in time and the conservation equation $\delta\sigma' = 0$ will be replaced by a dynamical one. On the other hand, if a solid is considered, an anisotropic stress is present and the two Bardeen potentials are not equal, so that

$$ \tau - \psi = 4M_{pl}^2\pi_L, \quad (8.16) $$

where $\pi_L$ is the longitudinal mode of the spatial Stückelberg perturbations, $\Phi^i = x^i + \partial_i\pi_L$, see also (7.20). We plan to do a complete analysis of both effects in a forthcoming paper.

## 1. Gravitational waves

The dynamics of gravitational waves is governed by (7.22). Without loss of generality, we can choose $N(t) = a(t)$, so that in Fourier space we have the action

$$ \frac{M_{pl}^2}{2} \int dt d^3 k a^2 [\chi'_{ij}\chi'_{ij} + (2M_{pl}^2 + k^2)\chi_{ij}\chi_{ij}]. \quad (8.17) $$

The resulting propagation equation is

$$ \chi_{ij}'' + 2H\chi_{ij}' + (k^2 - M_{pl}^2)\chi_{ij} = 0. \quad (8.18) $$

The difference with respect to the standard expression is the presence of the time-dependent effective mass squared, $M_{pl}^2$,

<footnote n="14">¹⁴Notice that in a more general case, where $c_s^2$ is not negligible, there is an interplay between $c_s^2$ and $c_b^2$ that makes the growth condition more complicated that just $c_b^2 < 1$.</footnote>
which is nonvanishing for solids and supersolids. As we had anticipated, this mass also controls the difference of the Bardeen potentials, see (8.16). We point out that, at leading order in derivatives, $M_2^2$ induces a modification of the speed of propagation (the dispersion relation) of gravitational waves, which becomes scale-dependent. Interestingly, it was shown in [16] that at next-to-leading order, the propagation speed changes in a scale-independent (but time-dependent) way in the case of VDiff invariant media made of just the three spatial Stückelberg fields. The same kind of modification is expected to occur at NLO for a generic medium, leading to a generic propagation equation of the form,

$$ \chi_{ij}'' + 2Q_1 \mathcal{H} \chi_{ij}' + (Q_2 k^2 - \tilde{M}_2^2) \chi_{ij} = 0, \quad (8.19) $$

where $\tilde{M}_2^2$, $Q_1$ and $Q_2$ are time-dependent functions and we expect $M_2^2/\tilde{M}_2^2 \sim 1 + O(H/\Lambda)$ and also $Q_1, Q_2 \sim 1 + O(H/\Lambda)$, where $\Lambda$ denotes the cutoff of the effective theory.

## IX. CONCLUSIONS AND OUTLOOK

The main result of this work is that a vast class of modified gravity models, including MG, can be interpreted as self-gravitating continuous media. The low-energy dynamics of these systems is described by the EFT of four scalar fields $\Phi^A$ that respect shift symmetries. Due to these symmetries, the fields are derivatively coupled between them and are minimally coupled to gravity at leading order in derivatives. At this order, the action is a functional of ten independent scalars encoded in the induced four-dimensional (inverse) metric $C^{AB} = g^{\mu\nu} \partial_\mu \Phi^A \partial_\nu \Phi^B$.

The four scalar fields can be interpreted as the (comoving) coordinates of the medium. They can be conveniently split into three spatial coordinates $\Phi^a$, $a = 1, 2, 3$, and a temporal coordinate $\Phi^0$. Using diffeomorphism invariance, the EFT can be examined in the unitary gauge, where the scalar fields are “frozen” to be coincident with the spacetime coordinates. With this choice of spacetime coordinates, the induced metric is $C^{AB} = g^{AB}$. This sends all the dynamics of the medium into the gravitational metric and allows us to make direct contact with the traditional framework of MG. The inverse path (going from MG to continuous media) can be travelled by using the well-known Stückelberg “trick,” which allows us to write in a diffeomorphism-invariant way any theory in which diffeomorphisms appear to be broken. Indeed, the four scalar fields $\Phi^A$ of our EFT can also be seen as the four Stückelberg fields of a diffeomorphism-invariant formulation of MG. We have discussed in Sec VI how to move from one picture to the other.

We have provided a comprehensive and systematic classification of massive gravity and modified gravity models in terms of symmetries and low-energy degrees

of freedom, establishing a direct connection to the theory of (self-gravitating) continuous media and the pull-back formalism. In doing so, we have determined also which are the propagating degrees of freedom for each case (as a function of the symmetries) and obtained the Lagrangians that are necessary for the study of (linear) cosmological perturbations at large scales, highlighting which operators affect the different kinds of observables.

The mechanical and thermodynamic properties of these continuous media (or modified/massive gravity) models depend crucially on extra symmetries that can be imposed on the scalar sector and may restrict drastically the number of allowed operators. The minimal assumption we have used in this work is that the action is invariant under (spatial) SO(3) rotations of the fields $\Phi^a$, $a = 1, 2, 3$. This reduces the number of allowed operators from ten to nine at leading order in derivatives; see Table I. Imposing further assumptions, such as symmetries that relate the spatial and temporal Stückelberg sectors—see Table II—the number of independent operators can be decreased even down to a single one. For instance, volume preserving internal spatial diffeomorphisms selects just three operators and the resulting action describes superfluids.

Another way to simplify the action is to make it depend (by assumption) only on either the spatial or temporal Stückelberg fields; see Table III. The models that are allowed in those cases can describe self-gravitating solids, perfect fluids and superfluids. In the general SO(3) spatially symmetric case, which contains the four scalars $\Phi^A$, the medium shares some of its properties with both superfluids and solids, and therefore is called a supersolid. Therefore, a general MG model is interpreted a supersolid propagating in spacetime.

Working in the unitary gauge and expanding a generic example of our EFT around Minkowski spacetime, it is straightforward to see that the resulting action for the metric fluctuations will not respect Lorentz symmetry. This is simply a natural consequence of the presence of the medium and the fact that we have assumed the operator content to be dictated by a broad symmetry: internal spatial SO(3) rotations. Nevertheless, it is possible to choose a combination of operators such that the resulting action for metric perturbations is Lorentz invariant, if one wishes to do so [13].

The correspondence between self-gravitating media and massive/modified gravity is intriguing and deserves further study. In the context of cosmology, it can be relevant for a deeper understanding of the properties of dark matter and dark energy. Moreover, the systematic EFT framework presented in this work can aid to devise a program for the interpretation of the data coming from future probes. From this perspective, the key idea is that the symmetry properties of a (coarse-grained) cosmological medium could be imprinted in the large-scale structure of the Universe, including the cosmic microwave background.
We plan to analyze the cosmology of these models in a future publication, by studying in further detail the signatures that would allow to distinguish these models from the standard ΛCDM model [79]. In the present work, we have initiated this program (in Secs. VII and VIII) by determining the number of propagating degrees of freedom in FLRW for various representative examples of continuous media (or models of massive/modified gravity), as summarized in Table V, and by studying analytically some representative examples. In particular, in Sec. VIII we have considered the modification to the growth function of matter in nonbarotropic perfect fluids and the modifications to the propagation of gravitational waves. Other avenues that are also worth exploring can be opened if different internal symmetries—such as supersymmetry [27] or Galileon symmetry [98]—are assumed instead of spatial SO(3).

## ACKNOWLEDGMENTS

We thank B. Bellazzini and A. Ramos for discussions. We thank D. Blas and S. Sibiryakov for useful comments on the draft. D. C. and L. P. would also like to thank F. Nesti for enlightening discussions during previous works. The work of G. B. is funded by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 656794. G. B. thanks the Deutsches Elektronen Synchrotron (DESY) theory group for hospitality while part of this work was developed and the German Science Foundation (DFG) for funding through the Collaborative Research Center (SFB) 676 "Particles, Strings and the Early Universe." G. B. thanks as well the CERN Theoretical Physics Department for hospitality while this work was completed. G. B. and L. P. would like to thank Nordita (Nordic Institute for Theoretical Physics) for hosting them during a fruitful workshop. L. P. thanks the Institute de Physique Théorique IPhT CEA-Saclay for hospitality during the preparation of this work.

## APPENDIX A: CONSERVED CURRENTS AND CHARGES

Each of the systems we consider in this paper is characterized by some conserved currents and their associated charges. These currents and charges have special significance for the thermodynamic and dynamical interpretation of the different media. Generically, the currents can be broadly classified in two different types. First, we have those currents which, by virtue of Noether's theorem, are due to the symmetries that define the Lagrangians of the different models. These currents are grouped in sets of infinite dimension because the defining symmetries are (infinite) subgroups of the internal diffeomorphisms $\Phi^A \to \Psi^A(\Phi^B)$. Then, we also have currents that are conserved irrespectively of the equations of motion. An example of the latter is

$$J^\mu = -bu^\mu, \qquad (A1)$$

where $u^\mu$ was defined in (3.4). This current is covariantly conserved off-shell and exists for all the models we consider (regardless of their symmetries) that involve the spatial Stückelberg fields $\Phi^a$. It has often been interpreted as the entropy current for fluid actions of the form $S = \int d^4xU(b)$, see e.g. [10,24], although this is not the only possible interpretation, see e.g. [99].

Independently of whether a current is conserved due to a symmetry of the action or for another reason, an associated charge that is conserved in time can be defined. In general, given a current $J^\mu$ that satisfies

$$\nabla_\mu J^\mu = 0, \qquad (A2)$$

the corresponding time-conserved charge is

$$Q = \int d^3x \sqrt{-g} J^0 = \int d^3\Phi \frac{J^0}{bu^0}, \qquad (A3)$$

which can be proven using the divergence theorem and assuming that the fields fall of quickly enough at infinity. For instance, in the particular case of the current (A1), the charge is

$$V_3 = - \int d^3\Phi \propto \epsilon_{ijk} \epsilon^{0a\beta\gamma} \int d^3x \partial_\alpha \Phi^i \partial_\beta \Phi^j \partial_\gamma \Phi^k \qquad (A4)$$

and, clearly, $\dot{V}_3 = 0$ due to the antisymmetric character of $\epsilon^{\mu a\beta\gamma}$. The physical interpretation of this result is that the flux lines of the medium are neither created nor destroyed, which is a topological statement. In other words, the conservation of $V_3$ means that these continuous media have neither flux sources nor sinks.

### 1. General volume currents

Let us first consider conserved currents that are independent of the internal symmetries. We just encountered one example of them in (A1). This type of current can be generalized in a simple way with a permutation of the fields $\Phi^A$. Concretely, the four currents

$$J_D^\mu = \frac{\epsilon^{\mu\alpha\beta\gamma}}{6\sqrt{-g}} \epsilon_{ABCD} \partial_\alpha \Phi^A \partial_\beta \Phi^B \partial_\gamma \Phi^C, \qquad D = 0, 1, 2, 3 \qquad (A5)$$

are covariantly conserved, and $J_0^\mu$ is precisely (A1). In fact, any four-vector of the form

$$\mathcal{J}_D^\mu [h(\Phi^A)] = h(\Phi^A) J_D^\mu, \qquad (A6)$$

where $h$ is an arbitrary function of the fields $\Phi^A$, satisfies
$$ \nabla_{\mu} \mathcal{J}_{D}^{\mu}[h(\Phi^{A})] = 4 \frac{\partial h(\Phi^{A})}{\partial \Phi^{D}} \det \left( \frac{\partial \Phi^{B}}{\partial x^{\mu}} \right). \quad (A7) $$

This allows to construct further conserved currents by simply multiplying any of the currents $J_D^\mu$ with a function $h$ that does not depend on the field $\Phi^D$. The relation (A7) can be easily proven using the identity $\sqrt{-g}\nabla_\mu U^\mu = \partial_\mu(\sqrt{-g}U^\mu)$, valid for any vector $U^\mu$, and the antisymmetry of $J_D^\mu$ under the permutation of two fields $\Phi^A$.

## 2. Noether currents

Let us now consider the currents coming from the different symmetries that we have considered in the paper. First of all, we have $SO(3)_S$ and the internal shifts (2.1). The first one of these two symmetries leads to the conservation of angular momentum in the internal spatial manifold. The conservation of the currents associated to the shift symmetries is nothing but the equations of motion of the Stückelber fields, see (2.5). We will now enumerate the conserved currents due to the additional symmetries of Sec. IVA.

### a. Four-dimensional media

(i) If the action is invariant under $V_S$Diff, the volume-preserving spatial diffeomorphisms in the space of Stückelberg fields, the operators $b$, $Y$ and $X$ are selected and the resulting LO medium is a super-fluid. The currents associated to this symmetry are known to be related to vorticity conservation, see for instance [9,10]. These currents are of the form

$$ J_{(\varepsilon)}^{\mu} = -b U_b (\mathbf{B}^{-1})^{cd} \partial^{\mu} \Phi^d \varepsilon^c (\Phi^j), \quad (A8) $$

where $\varepsilon^a(\Phi^j)$ satisfies $\partial\varepsilon^a/\partial\Phi^a = 0$. A basis of vorticity charges (in which any other vorticity charge can be expressed) is constructed choosing $\varepsilon^a = \varepsilon_{abc}\alpha_b \partial^3(\Phi - \tilde{\Phi})/\partial\tilde{\Phi}^c$, with constant $\alpha_b$, see e.g. [10].

(ii) The covariantly conserved Noether currents associated to the symmetry $\Phi^a \to \Phi^a + f^a(\Phi^0)$, which selects the operators $X$ and $w_n$, are

$$ J_f^\mu = \sum_{a=1}^{3} f^a(\Phi^0) \sum_{n=1}^{3} n U_{w_n} [W^{n-1} \cdot j]_a^\mu, \quad (A9) $$

where we define

$$ j_a^\mu = (\partial^\mu \Phi^a - C^a X^{-1} \partial^\mu \Phi^0). \quad (A10) $$

Each set of three functions $\{f^1, f^2, f^3\}$ of the temporal Stückelberg field defines a different current. By using the unitary gauge expressions (6.5) we interpret $\Phi^a \to \Phi^a + f^a(\Phi^0)$ as a transformation

that does not change the slicing of spacetime but modifies the threading with a time-dependent shift: $N \to N$, $N^i \to N^i + N^{-1}\partial f^i / \partial \Phi^0$. Notice the spatial part of the metric, $\gamma_{ij}$, is invariant under $\Phi^a \to \Phi^a + f^a(\Phi^0)$.

(iii) In the case of the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^a)$, the symmetric LO scalar operators are $Y$ and $\tau_n$, and the covariantly conserved currents are

$$ J_f^\mu = f(\Phi^a) Y U_Y u^\mu. \quad (A11) $$

It seems natural to interpret this set of currents—recall that there is a current for each choice of the function $f$—as transporting charges of the fluid in the direction of the four-velocity $u^\mu$. Notice that these currents are parallel to the current (A1). Therefore, if (A1) is interpreted as the entropy current, the currents (A11) carry charges that flow with the entropy.

(iv) If the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^0)$ is imposed, the LO action depends on the operators $\tau_n$ and $O_{\alpha\beta n}$. The operators $\tau_n$ do not contribute to the conserved currents, since they do not contain $\Phi^0$. The currents in this case are

$$ J_f^\mu = \frac{f(\Phi^0)}{Y^2} \sum_{\alpha, \beta, n} U_{\alpha\beta n} \quad (A12) $$

where the sums extend over all the $O_{\alpha\beta n}$ operators and

$$ \begin{aligned} J_{\alpha\beta\gamma}^{\mu} = & \alpha O_{(\alpha-1)\beta n} \partial^{\mu} \phi^{0} - (\alpha + \beta) Y O_{\alpha\beta n} u^{\mu} \\ & + \beta O_{\alpha(\beta-1)n} (\mathbf{B}^{n})^{ab} C^{0a} \partial^{\mu} \Phi^{b}. \end{aligned} \quad (A13) $$

As we discussed in Sec. IVA, the symmetry $\Phi^0 \to \Phi^0 + f(\Phi^0)$ allows at LO an infinite number of different scalars $O_{\alpha\beta n}$. This is because the exponents $\alpha$ and $\beta$ appearing in the definition of $O_{\alpha\beta n}$ are arbitrary real numbers (whereas $n$ can take the values 1, 2 or 3). In consequence the summatory $\sum_{\alpha, \beta, n}$ has to be understood to extend over all the possible values of these parameters, unless some restriction is enforced on them.

### b. Media with reduced internal dimensionality

(i) We can also consider the case in which the Lagrangian depends only on the temporal Stückelberg field, although, as we explained is Sec. VA 3, we do not know of any symmetry that forbids the spatial Stückelberg fields at LO (and allows them at higher orders). In this case the master function $U$ depends exclusively on $X$ at LO and the relevant symmetry is the invariance under $\Phi^0 \to \Phi^0 + c^0$, with $\partial_\mu c^0 = 0$. This simply tells us that the current
$$ \mathcal{X}^{\mu} = \sqrt{-X} U_{\chi} V^{\mu} \qquad (A14) $$

is covariantly conserved. The conservation of $\mathcal{X}^{\mu}$ is precisely the equation of motion for $\Phi^0$, given the master function $U(X)$.

(ii) Analogously, we can consider the case in which only the spatial Stückelberg fields are present. Since we are assuming always an $SO(3)_s$ symmetry this corresponds to the standard solid $U(\tau_n)$. Notice, again, that we know of no symmetry that prevents the appearance of $\Phi^0$ at LO but reintroduces this field at higher orders. Aside from the conservation of internal angular momentum, due to $SO(3)_s$, there are also the currents associated to the internal translational invariance $\Phi^a \to \Phi^a + c^a$, with $\partial_\mu c^a = 0$. This symmetry generates in this case the set of covariantly conserved currents

$$ t^{\mu} = c^{a} \sum_{n=1}^{3} n U_{\tau_{n}} (\mathbf{B}^{n-1})^{ab} \partial^{\mu} \Phi^{b}, \qquad (A15) $$

a base of which is obtained by setting, alternatively, $c^a$ to 0 or 1 for $a = 1, 2$ or $3$.

The currents associated to $SO(3)_s$ invariance have a similar structure. In fact, they can be formally obtained from the previous ones simply replacing $c^a$ by $r^{ab}\Phi^b$, where $r^{ab}$ belongs to the Lie algebra of $SO(3)_s$, and so is antisymmetric. These are just the generators of the angular momentum in the internal spatial manifold. Therefore, the currents

$$ J_r^\mu = r^{ab} \Phi^b \sum_{n=1}^{3} n U_{\tau_n} (\mathbf{B}^{n-1})^{ab} \partial^\mu \Phi^b, \qquad (A16) $$

are also conserved in this case. A basis for this currents can be obtained choosing the matrices $r^{ab}$ as the standard $L_x, L_y, L_z$ generators. The conservation of the currents of this basis are the equations of motion for the spatial Stückelberg fields.

## APPENDIX B: A SHORTCUT TO COUNT DEGREES OF FREEDOM

We can determine the number of d.o.f. by studying directly to the equations of motion of the Stückelberg fields (2.5). From the action (2.3) we know that there are at most

six d.o.f., that is, two from the metric plus those coming from the scalar fields ($\le 4$). Equation (2.5) can be rewritten in the following form (using the notation $U_{AB} = \partial U/\partial C^{AB}$):

$$ \begin{aligned} \mathcal{L}_A &= [U_{AB}g^{\mu\nu} + (U_{AD,CB} + U_{AD,BC})\nabla^{\mu}\Phi^C\nabla^{\nu}\Phi^D]\nabla_{\mu}\nabla_{\nu}\Phi^B \\ &\equiv \mathcal{K}_{AB}^{\mu\nu} \nabla_{\mu}\nabla_{\nu}\Phi^B = 0, \end{aligned} \qquad (B1) $$

where $\mathcal{K}_{AB}^{\mu\nu}$ contains only first derivatives of the fields. From this expression we see that the second-order time derivatives of $\Phi^A$ are always proportional to the matrix kernel $\mathcal{K}_{AB}^{00}$, thus the number of Stückelberg propagating d.o.f. is

$$ \Delta = \mathrm{Rank}[\mathcal{K}_{AB}^{00}]. \qquad (B2) $$

To simplify the analysis it is convenient to use the phonons (2.8) and study the system perturbatively. The coefficient of second-order time derivatives of the perturbations are determined by $\bar{\mathcal{K}}_{AB}^{00} \equiv \mathcal{K}_{AB}^{00}|_{\Phi^A=x^A}$. From (B1) we have that $(\bar{\mathcal{K}}_{AB}^{00} + \mathcal{O}(\partial\pi))\hat{x}^B = 0$ where

$$ \bar{\mathcal{K}}_{AB}^{00} = \bar{L}_{AB}g^{00} + (\bar{L}_{AD,CB} + \bar{L}_{AD,BC})g^{0C}g^{0D}, \qquad (B3) $$

where the overbars mean that the corresponding quantities are evaluated on $\Phi^A = x^A$. It is easy to see that $\bar{\mathcal{K}}_{AB}^{00}$ gives the number of d.o.f. also at the nonperturbative level. Indeed, by going to the unitary gauge we have that Rank[$\mathcal{K}_{AB}^{00}$] = $\Delta$ and so the counting of d.o.f. is background independent. For $SO(3)_s$ invariant Lagrangians one finds the same results obtained applying the Hamiltonian formalism [11–13], that is,

(i) for $U = L(g^{AB})$ we get four scalar d.o.f.,

(ii) for $U = L(g^{ab}) + \sqrt{-g^{00}}\tilde{L}(\gamma^{ab})$ we get three scalar d.o.f.,

(iii) for $L_1 = L(g^{00}, \gamma^{ab})$ we get one scalar d.o.f.,

(iv) for $L_0 = \sqrt{-g^{00}}\tilde{L}(\gamma^{ab})$ we get zero scalar d.o.f.

The Hamiltonian analysis applied to field theories allows for a noninteger number of d.o.f. (see [11,13]) while the prescription of (B2) seems to evade such a problem. For example the naive Hamiltonian counting of d.o.f. for a master function of the form $U = \sum_a C^{0a}C^{0a}/C^{00}$ gives 3 + 1/2 d.o.f. in the Hamiltonian formalism and three d.o.f. when computed following (B2). This mismatch is being studied.
[1] A. G. Riess et al. (Supernova Search Team Collaboration), Astron. J. **116**, 1009 (1998); S. Perlmutter et al. (Supernova Cosmology Project Collaboration), Astrophys. J. **517**, 565 (1999).

[2] P. A. R. Ade et al. (Planck Collaboration), Astron. Astrophys. **594**, A14 (2016).

[3] P. A. R. Ade et al. (Planck Collaboration), Astron. Astrophys. **594**, A13 (2016).

[4] R. Laureijs et al. (Euclid Collaboration), arXiv:1110.3193; J. Green et al. (WFIRST), arXiv:1208.4012; M. Levi et al. (DESI Collaboration), arXiv:1308.0847; A. Abate et al. (LSST Dark Energy Science Collaboration), arXiv: 1211.0310.

[5] L. Amendola and S. Tsujikawa, *Dark Energy, Theory and Observations* (Cambridge University Press, Cambridge, England, 2010); E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D **15**, 1753 (2006); T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. **513**, 1 (2012).

[6] S. Weinberg, Rev. Mod. Phys. **61**, 1 (1989).

[7] S. L. Dubovsky, J. High Energy Phys. 10 (2004) 076.

[8] V. A. Rubakov and P. G. Tinyakov, Phys. Usp. **51**, 759 (2008).

[9] S. Dubovsky, T. Gregoire, A. Nicolis, and R. Rattazzi, J. High Energy Phys. 03 (2006) 025.

[10] G. Ballesteros and B. Bellazzini, J. Cosmol. Astropart. Phys. 04 (2013) 001.

[11] D. Comelli, M. Crisostomi, F. Nesti, and L. Pilo, Phys. Rev. D **86**, 101502 (2012).

[12] D. Comelli, F. Nesti, and L. Pilo, J. High Energy Phys. 07 (2013) 161.

[13] D. Comelli, F. Nesti, and L. Pilo, J. Cosmol. Astropart. Phys. 11 (2014) 018.

[14] G. Gubitosi, F. Piazza, and F. Vernizzi, J. Cosmol. Astropart. Phys. 02 (2013) 032.

[15] J. Bhattacharya, S. Bhattacharyya, and M. Rangamani, J. High Energy Phys. 02 (2013) 153.

[16] G. Ballesteros, J. Cosmol. Astropart. Phys. 03 (2015) 001.

[17] R. M. Wald, *General Relativity* (University of Chicago Press, Chicago, 1984).

[18] E. C. G. Stückelberg, Helv. Phys. Acta **11**, 225 (1938); **11**, 299 (1938); **11**, 312 (1938).

[19] H. Ruegg and M. Ruiz-Altaba, Int. J. Mod. Phys. A **19**, 3265 (2004).

[20] M. Fierz and W. Pauli, Proc. R. Soc. A **173**, 211 (1939).

[21] H. Leutwyler, Helv. Phys. Acta **70**, 275 (1997).

[22] H. Leutwyler, Phys. Rev. D **49**, 3033 (1994).

[23] S. Endlich, A. Nicolis, R. Rattazzi, and J. Wang, J. High Energy Phys. 04 (2011) 102.

[24] S. Dubovsky, L. Hui, A. Nicolis, and D. T. Son, Phys. Rev. D **85**, 085029 (2012).

[25] S. Dubovsky, L. Hui, and A. Nicolis, Phys. Rev. D **89**, 045016 (2014).

[26] A. Nicolis, arXiv:1108.2513.

[27] C. Hoyos, B. Keren-Zur, and Y. Oz, J. High Energy Phys. 11 (2012) 152.

[28] S. Endlich, A. Nicolis, and J. Wang, J. Cosmol. Astropart. Phys. 10 (2013) 011.

[29] S. Endlich, A. Nicolis, R. A. Porto, and J. Wang, Phys. Rev. D **88**, 105001 (2013).

[30] A. Nicolis, R. Penco, and R. A. Rosen, Phys. Rev. D **89**, 045002 (2014).

[31] G. Ballesteros, B. Bellazzini, and L. Mercolli, J. Cosmol. Astropart. Phys. 05 (2014) 007.

[32] B. Gripaios and D. Sutherland, Phys. Rev. Lett. **114**, 071601 (2015).

[33] L. V. Delacrétaz, A. Nicolis, R. Penco, and R. A. Rosen, Phys. Rev. Lett. **114**, 091601 (2015).

[34] B. Carter and H. Quintana, Proc. R. Soc. A **331**, 57 (1972).

[35] B. Carter, *Relativistic Fluid Dynamics* (Springer, New York, 1989), pp. 1–64.

[36] N. Andersson and G. L. Comer, Living Rev. Relativ. **10**, 1 (2007).

[37] J. D. Brown, Classical Quantum Gravity **10**, 1579 (1993).

[38] G. L. Comer and D. Langlois, Classical Quantum Gravity **10**, 2317 (1993).

[39] G. L. Comer and D. Langlois, Classical Quantum Gravity **11**, 709 (1994).

[40] A. Gruzinov, Phys. Rev. D **70**, 063518 (2004).

[41] D. Blas, M. M. Ivanov, and S. Sibiryakov, J. Cosmol. Astropart. Phys. 10 (2012) 057.

[42] C. Skordis, A. Pourtsidou, and E. J. Copeland, Phys. Rev. D **91**, 083537 (2015).

[43] L. Berezhiani and J. Khoury, Phys. Rev. D **92**, 103510 (2015).

[44] N. Bartolo, S. Matarrese, M. Peloso, and A. Ricciardone, J. Cosmol. Astropart. Phys. 08 (2013) 022.

[45] S. Endlich, B. Horn, A. Nicolis, and J. Wang, Phys. Rev. D **90**, 063506 (2014).

[46] D. T. Son, Phys. Rev. Lett. **94**, 175301 (2005).

[47] B. Carter and H. Quintana, Proc. R. Soc. A **331**, 57 (1972); J. Kijowski and G. Magli, J. Geom. Phys. **9**, 207 (1992); J. Kijowski and G. Magli, arXiv:hep-th/9411212.

[48] L. D. Landau and E. M. Lifshitz, *Elasticity Theory: Course of Theoretical Physics* (Pergamon Press, Oxford, 1970), Vol. 7.

[49] A. F. Bower, *Applied Mechanics of Solids* (CRC Press, Boca Raton, 2009).

[50] B. Carter, Commun. Math. Phys. **30**, 261 (1973).

[51] M. Bucher and D. N. Spergel, Phys. Rev. D **60**, 043505 (1999).

[52] G. Ballesteros, D. Comelli, and L. Pilo, Phys. Rev. D **94**, 025034 (2016).

[53] S. Matarrese, Proc. R. Soc. A **401**, 53 (1985).

[54] T. Jacobson and D. Mattingly, Phys. Rev. D **64**, 024028 (2001).

[55] D. Blas, O. Pujolas, and S. Sibiryakov, J. High Energy Phys. 04 (2011) 018.

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

[57] C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, J. High Energy Phys. 03 (2008) 014.

[58] E. Kim and M. H. W. Chan, Nature (London) **427**, 225 (2004).

[59] M. Boninsegni and N. V. Prokofiev, Rev. Mod. Phys. **84**, 759 (2012).

[60] C. G. Callan, Jr., S. R. Coleman, J. Wess, and B. Zumino, Phys. Rev. **177**, 2247 (1969).

[61] D. T. Son, arXiv:hep-ph/0204199.

[62] D. G. Boulware and S. Deser, Phys. Lett. **40B**, 227 (1972).
[63] C. de Rham, G. Gabadadze, and A. J. Tolley, Phys. Rev. Lett. **106**, 231101 (2011).

[64] S. F. Hassan and R. A. Rosen, Phys. Rev. Lett. **108**, 041101 (2012).

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

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

[67] C. de Rham, A. J. Tolley, and S. Y. Zhou, J. High Energy Phys. 04 (2016) 188.

[68] C. de Rham, A. J. Tolley, and S. Y. Zhou, Phys. Lett. B **760**, 579 (2016).

[69] H. van Dam and M. J. G. Veltman, Nucl. Phys. **B22**, 397 (1970); Y. Iwasaki, Phys. Rev. D **2**, 2255 (1970); V. I. Zakharov, ZhETF Pis. Red. **12**, 447 (1970) [JETP Lett. **12**, 198 (1971)].

[70] A. I. Vainshtein, Phys. Lett. **39B**, 393 (1972).

[71] E. Babichev, C. Deffayet, and R. Ziour, Phys. Rev. Lett. **103**, 201102 (2009); N. Kaloper, A. Padilla, and N. Tanahashi, J. High Energy Phys. 10 (2011) 148; G. Chkareuli and D. Pirtskhalava, Phys. Lett. B **713**, 99 (2012).

[72] D. Blas and E. Lim, Int. J. Mod. Phys. D **23**, 1443009 (2014).

[73] V. A. Rubakov, arXiv:hep-th/0407104.

[74] D. Comelli, F. Nesti, and L. Pilo, J. Cosmol. Astropart. Phys. 05 (2014) 036.

[75] W. Siegel, Phys. Rev. D **49**, 4144 (1994).

[76] K. Hinterbichler, Rev. Mod. Phys. **84**, 671 (2012).

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

[78] D. Blas and S. Sibiryakov, Zh. Eksp. Teor. Fiz. **147**, 578 (2015) [J. Exp. Theor. Phys. **120**, 509 (2015)].

[79] G. Ballesteros, D. Comelli, and L. Pilo (to be published).

[80] L. Amendola et al., arXiv:1606.00180.

[81] L. M. Wang and P. J. Steinhardt, Astrophys. J. **508**, 483 (1998).

[82] E. V. Linder, Phys. Rev. D **72**, 043529 (2005).

[83] E. Bertschinger and P. Zukin, Phys. Rev. D **78**, 024015 (2008).

[84] G. Ballesteros and A. Riotto, Phys. Lett. B **668**, 171 (2008).

[85] L. Perivolaropoulos, arXiv:0811.4684.

[86] R. Bean and O. Dore, Phys. Rev. D **69**, 083503 (2004).

[87] S. Hannestad, Phys. Rev. D **71**, 103519 (2005).

[88] M. Takada, Phys. Rev. D **74**, 043505 (2006).

[89] R. de Putter, D. Huterer, and E. V. Linder, Phys. Rev. D **81**, 103513 (2010).

[90] G. Ballesteros and J. Lesgourgues, J. Cosmol. Astropart. Phys. 10 (2010) 014.

[91] D. Sapone, Int. J. Mod. Phys. A **25**, 5253 (2010).

[92] J. M. Bardeen, Phys. Rev. D **22**, 1882 (1980).

[93] L. Amendola, M. Kunz, and D. Sapone, J. Cosmol. Astropart. Phys. 04 (2008) 013.

[94] S. F. Daniel, R. R. Caldwell, A. Cooray, and A. Melchiorri, Phys. Rev. D **77**, 103513 (2008).

[95] L. Amendola, G. Ballesteros, and V. Pettorino, Phys. Rev. D **90**, 043009 (2014).

[96] M. Raveri, C. Baccigalupi, A. Silvestri, and S. Y. Zhou, Phys. Rev. D **91**, 061501 (2015).

[97] I. D. Saltas, I. Sawicki, L. Amendola, and M. Kunz, Phys. Rev. Lett. **113**, 191101 (2014).

[98] A. Nicolis, R. Penco, F. Piazza, and R. Rattazzi, J. High Energy Phys. 06 (2015) 155.

[99] B. F. Schutz, Phys. Rev. D **2**, 2762 (1970).